Actual source code: ex28.c
3: static char help[] = "Solves 1D wave equation using multigrid.\n\n";
5: #include petscda.h
6: #include petscksp.h
7: #include petscdmmg.h
15: int main(int argc,char **argv)
16: {
18: PetscInt i;
19: DMMG *dmmg;
20: PetscReal norm;
21: DA da;
23: PetscInitialize(&argc,&argv,(char *)0,help);
25: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
26: DACreate1d(PETSC_COMM_WORLD,DA_XPERIODIC,-3,2,1,0,&da);
27: DMMGSetDM(dmmg,(DM)da);
28: DADestroy(da);
30: DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
32: ComputeInitialSolution(dmmg);
34: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
35: for (i=0; i<1000; i++) {
36: DMMGSolve(dmmg);
37: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
38: }
39: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
40: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
41: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
42: /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm); */
44: DMMGDestroy(dmmg);
45: PetscFinalize();
47: return 0;
48: }
52: PetscErrorCode ComputeInitialSolution(DMMG *dmmg)
53: {
55: PetscInt mx,col[2],xs,xm,i;
56: PetscScalar Hx,val[2];
57: Vec x = DMMGGetx(dmmg);
60: DAGetInfo(DMMGGetDA(dmmg),0,&mx,0,0,0,0,0,0,0,0,0);
61: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
62: DAGetCorners(DMMGGetDA(dmmg),&xs,0,0,&xm,0,0);
63:
64: for(i=xs; i<xs+xm; i++){
65: col[0] = 2*i; col[1] = 2*i + 1;
66: val[0] = val[1] = PetscSinScalar(((PetscScalar)i)*Hx);
67: VecSetValues(x,2,col,val,INSERT_VALUES);
68: }
69: VecAssemblyBegin(x);
70: VecAssemblyEnd(x);
71: return(0);
72: }
73:
76: PetscErrorCode ComputeRHS(DMMG dmmg,Vec b)
77: {
79: PetscInt mx;
80: PetscScalar h;
83: DAGetInfo((DA)dmmg->dm,0,&mx,0,0,0,0,0,0,0,0,0);
84: h = 2.0*PETSC_PI/((mx));
85: VecCopy(dmmg->x,b);
86: VecScale(b,h);
87: return(0);
88: }
92: PetscErrorCode ComputeMatrix(DMMG dmmg,Mat J,Mat jac)
93: {
94: DA da = (DA)dmmg->dm;
96: PetscInt i,mx,xm,xs;
97: PetscScalar v[7],Hx;
98: MatStencil row,col[7];
99: PetscScalar lambda;
101: PetscMemzero(col,7*sizeof(MatStencil));
102: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
103: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
104: DAGetCorners(da,&xs,0,0,&xm,0,0);
105: lambda = 2.0*Hx;
106: for(i=xs; i<xs+xm; i++){
107: row.i = i; row.j = 0; row.k = 0; row.c = 0;
108: v[0] = Hx; col[0].i = i; col[0].c = 0;
109: v[1] = lambda; col[1].i = i-1; col[1].c = 1;
110: v[2] = -lambda;col[2].i = i+1; col[2].c = 1;
111: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
113: row.i = i; row.j = 0; row.k = 0; row.c = 1;
114: v[0] = lambda; col[0].i = i-1; col[0].c = 0;
115: v[1] = Hx; col[1].i = i; col[1].c = 1;
116: v[2] = -lambda;col[2].i = i+1; col[2].c = 0;
117: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
118: }
119: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
120: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
121: MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
122: return 0;
123: }