Actual source code: ex34.c

```  1: /*T
2:    Concepts: KSP^solving a system of linear equations
3:    Concepts: KSP^Laplacian, 3d
4:    Processors: n
5: T*/

7: /*
8: Laplacian in 3D. Modeled by the partial differential equation

10:    div  grad u = f,  0 < x,y,z < 1,

12: with pure Neumann boundary conditions

14:    u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.

16: The functions are cell-centered

18: This uses multigrid to solve the linear system

20:        Contributed by Jianming Yang <jianming-yang@uiowa.edu>
21: */

23: static char help[] = "Solves 3D Laplacian using multigrid.\n\n";

25:  #include petscda.h
26:  #include petscksp.h
27:  #include petscmg.h
28:  #include petscdmmg.h

33: typedef enum {DIRICHLET, NEUMANN} BCType;

35: typedef struct {
36:   BCType        bcType;
37: } UserContext;

41: int main(int argc,char **argv)
42: {
43:   DMMG           *dmmg;
44:   DA             da;
45:   UserContext    user;
46:   PetscReal      norm;
47:   const char     *bcTypes[2] = {"dirichlet","neumann"};
49:   PetscInt       l,bc;

51:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
52:   PetscScalar    Hx,Hy,Hz;
53:   PetscScalar    ***array;

56:   PetscInitialize(&argc,&argv,(char *)0,help);
57:
58:   DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
59:   DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
60:   DASetInterpolationType(da, DA_Q0);

62:   DMMGSetDM(dmmg,(DM)da);
64:   for (l = 0; l < DMMGGetLevels(dmmg); l++) {
65:     DMMGSetUser(dmmg,l,&user);
66:   }
67:
68:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMMG");
69:   bc          = (PetscInt)NEUMANN;
70:   PetscOptionsEList("-bc_type","Type of boundary condition","ex34.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
71:   user.bcType = (BCType)bc;
72:   PetscOptionsEnd();
73:
74:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
75:   if (user.bcType == NEUMANN) {
76:     DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
77:   }

79:   DMMGSolve(dmmg);
80:
81:   MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
82:   VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
83:   VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
84:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
85:
86:   DAGetInfo(DMMGGetDA(dmmg), 0, &mx, &my, &mz, 0,0,0,0,0,0,0);
87:   Hx   = 1.0 / (PetscReal)(mx);
88:   Hy   = 1.0 / (PetscReal)(my);
89:   Hz   = 1.0 / (PetscReal)(mz);
90:   DAGetCorners(DMMGGetDA(dmmg),&xs,&ys,&zs,&xm,&ym,&zm);
91:   DAVecGetArray(DMMGGetDA(dmmg), DMMGGetx(dmmg), &array);

93:   for (k=zs; k<zs+zm; k++){
94:     for (j=ys; j<ys+ym; j++){
95:       for(i=xs; i<xs+xm; i++){
96:         array[k][j][i] -=
97:           PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
98:           PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
99:           PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
100:       }
101:     }
102:   }
103:   DAVecRestoreArray(DMMGGetDA(dmmg), DMMGGetx(dmmg), &array);
104:   VecAssemblyBegin(DMMGGetx(dmmg));
105:   VecAssemblyEnd(DMMGGetx(dmmg));

107:   VecNorm(DMMGGetx(dmmg),NORM_INFINITY,&norm);
108:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm);
109:   VecNorm(DMMGGetx(dmmg),NORM_1,&norm);
110:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
111:   VecNorm(DMMGGetx(dmmg),NORM_2,&norm);
112:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));

114:   DMMGDestroy(dmmg);
115:   PetscFinalize();
116:   return 0;
117: }

121: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
122: {
123:   DA             da = (DA)dmmg->dm;
124:   UserContext    *user = (UserContext *) dmmg->user;
126:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
127:   PetscScalar    Hx,Hy,Hz;
128:   PetscScalar    ***array;

131:   DAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0);
132:   Hx   = 1.0 / (PetscReal)(mx);
133:   Hy   = 1.0 / (PetscReal)(my);
134:   Hz   = 1.0 / (PetscReal)(mz);
135:   DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
136:   DAVecGetArray(da, b, &array);
137:   for (k=zs; k<zs+zm; k++){
138:     for (j=ys; j<ys+ym; j++){
139:       for(i=xs; i<xs+xm; i++){
140:         array[k][j][i] = 12*PETSC_PI*PETSC_PI
141:           *PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
142:           *PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
143:           *PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
144:           *Hx*Hy*Hz;
145:       }
146:     }
147:   }
148:   DAVecRestoreArray(da, b, &array);
149:   VecAssemblyBegin(b);
150:   VecAssemblyEnd(b);

152:   /* force right hand side to be consistent for singular matrix */
153:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
154:   if (user->bcType == NEUMANN)
155:     {
156:       MatNullSpace nullspace;
157:
158:       KSPGetNullSpace(dmmg->ksp,&nullspace);
159:       MatNullSpaceRemove(nullspace,b,PETSC_NULL);
160:     }

162:   return(0);
163: }

165:
168: PetscErrorCode ComputeJacobian(DMMG dmmg, Mat J,Mat jac)
169: {
170:   DA             da = (DA) dmmg->dm;
171:   UserContext    *user = (UserContext *) dmmg->user;
173:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
174:   PetscScalar    v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
175:   MatStencil     row, col[7];

178:   DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
179:   Hx    = 1.0 / (PetscReal)(mx);
180:   Hy    = 1.0 / (PetscReal)(my);
181:   Hz    = 1.0 / (PetscReal)(mz);
182:   HyHzdHx = Hy*Hz/Hx;
183:   HxHzdHy = Hx*Hz/Hy;
184:   HxHydHz = Hx*Hy/Hz;
185:   DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
186:   for (k=zs; k<zs+zm; k++)
187:     {
188:       for (j=ys; j<ys+ym; j++)
189:         {
190:           for(i=xs; i<xs+xm; i++)
191:             {
192:               row.i = i; row.j = j; row.k = k;
193:               if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1)
194:                 {
195:                   if (user->bcType == DIRICHLET)
196:                     {
197:                       SETERRQ(PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
198:                       v[0] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz);
199:                       MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
200:
201:                     }
202:                   else if (user->bcType == NEUMANN)
203:                     {
204:                       num = 0; numi=0; numj=0; numk=0;
205:                       if (k!=0)
206:                         {
207:                           v[num] = -HxHydHz;
208:                           col[num].i = i;
209:                           col[num].j = j;
210:                           col[num].k = k-1;
211:                           num++; numk++;
212:                         }
213:                       if (j!=0)
214:                         {
215:                           v[num] = -HxHzdHy;
216:                           col[num].i = i;
217:                           col[num].j = j-1;
218:                           col[num].k = k;
219:                           num++; numj++;
220:                         }
221:                       if (i!=0)
222:                         {
223:                           v[num] = -HyHzdHx;
224:                           col[num].i = i-1;
225:                           col[num].j = j;
226:                           col[num].k = k;
227:                           num++; numi++;
228:                         }
229:                       if (i!=mx-1)
230:                         {
231:                           v[num] = -HyHzdHx;
232:                           col[num].i = i+1;
233:                           col[num].j = j;
234:                           col[num].k = k;
235:                           num++; numi++;
236:                         }
237:                       if (j!=my-1)
238:                         {
239:                           v[num] = -HxHzdHy;
240:                           col[num].i = i;
241:                           col[num].j = j+1;
242:                           col[num].k = k;
243:                           num++; numj++;
244:                         }
245:                       if (k!=mz-1)
246:                         {
247:                           v[num] = -HxHydHz;
248:                           col[num].i = i;
249:                           col[num].j = j;
250:                           col[num].k = k+1;
251:                           num++; numk++;
252:                         }
253:                       v[num]   = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
254:                       col[num].i = i;   col[num].j = j;   col[num].k = k;
255:                       num++;
256:                       MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
257:                     }
258:                 }
259:               else
260:                 {
261:                   v[0] = -HxHydHz;                          col[0].i = i;   col[0].j = j;   col[0].k = k-1;
262:                   v[1] = -HxHzdHy;                          col[1].i = i;   col[1].j = j-1; col[1].k = k;
263:                   v[2] = -HyHzdHx;                          col[2].i = i-1; col[2].j = j;   col[2].k = k;
264:                   v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i;   col[3].j = j;   col[3].k = k;
265:                   v[4] = -HyHzdHx;                          col[4].i = i+1; col[4].j = j;   col[4].k = k;
266:                   v[5] = -HxHzdHy;                          col[5].i = i;   col[5].j = j+1; col[5].k = k;
267:                   v[6] = -HxHydHz;                          col[6].i = i;   col[6].j = j;   col[6].k = k+1;
268:                   MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
269:                 }
270:             }
271:         }
272:     }
273:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
274:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

276:   return(0);
277: }

```