Actual source code: ex34.c


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

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

  7: with pure Neumann boundary conditions

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

 11: The functions are cell-centered

 13: This uses multigrid to solve the linear system

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

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

 20: #include <petscdm.h>
 21: #include <petscdmda.h>
 22: #include <petscksp.h>

 24: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
 25: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);

 27: int main(int argc,char **argv)
 28: {
 29:   KSP            ksp;
 30:   DM             da;
 31:   PetscReal      norm;
 32:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,d,dof;
 33:   PetscScalar    Hx,Hy,Hz;
 34:   PetscScalar    ****array;
 35:   Vec            x,b,r;
 36:   Mat            J;

 38:   PetscInitialize(&argc,&argv,(char*)0,help);
 39:   dof  = 1;
 40:   PetscOptionsGetInt(NULL,NULL,"-da_dof",&dof,NULL);
 41:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 42:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,12,12,12,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,1,0,0,0,&da);
 43:   DMSetFromOptions(da);
 44:   DMSetUp(da);
 45:   DMDASetInterpolationType(da, DMDA_Q0);

 47:   KSPSetDM(ksp,da);

 49:   KSPSetComputeRHS(ksp,ComputeRHS,NULL);
 50:   KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
 51:   KSPSetFromOptions(ksp);
 52:   KSPSolve(ksp,NULL,NULL);
 53:   KSPGetSolution(ksp,&x);
 54:   KSPGetRhs(ksp,&b);
 55:   KSPGetOperators(ksp,NULL,&J);
 56:   VecDuplicate(b,&r);

 58:   MatMult(J,x,r);
 59:   VecAXPY(r,-1.0,b);
 60:   VecNorm(r,NORM_2,&norm);
 61:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

 63:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
 64:   Hx   = 1.0 / (PetscReal)(mx);
 65:   Hy   = 1.0 / (PetscReal)(my);
 66:   Hz   = 1.0 / (PetscReal)(mz);
 67:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
 68:   DMDAVecGetArrayDOF(da, x, &array);

 70:   for (k=zs; k<zs+zm; k++) {
 71:     for (j=ys; j<ys+ym; j++) {
 72:       for (i=xs; i<xs+xm; i++) {
 73:         for (d=0; d<dof; d++) {
 74:           array[k][j][i][d] -=
 75:             PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
 76:             PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
 77:             PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
 78:         }
 79:       }
 80:     }
 81:   }
 82:   DMDAVecRestoreArrayDOF(da, x, &array);
 83:   VecAssemblyBegin(x);
 84:   VecAssemblyEnd(x);

 86:   VecNorm(x,NORM_INFINITY,&norm);
 87:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)norm);
 88:   VecNorm(x,NORM_1,&norm);
 89:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
 90:   VecNorm(x,NORM_2,&norm);
 91:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
 92:   VecSum(x,&norm);
 93:   if (norm > 10000*PETSC_MACHINE_EPSILON) {
 94:     PetscPrintf(PETSC_COMM_WORLD,"Vector sum %g\n",(double)norm);
 95:   }

 97:   VecDestroy(&r);
 98:   KSPDestroy(&ksp);
 99:   DMDestroy(&da);
100:   PetscFinalize();
101:   return 0;
102: }

104: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
105: {
106:   PetscInt       d,dof,i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
107:   PetscScalar    Hx,Hy,Hz;
108:   PetscScalar    ****array;
109:   DM             da;
110:   MatNullSpace   nullspace;

113:   KSPGetDM(ksp,&da);
114:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,&dof,0,0,0,0,0);
115:   Hx   = 1.0 / (PetscReal)(mx);
116:   Hy   = 1.0 / (PetscReal)(my);
117:   Hz   = 1.0 / (PetscReal)(mz);
118:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
119:   DMDAVecGetArrayDOFWrite(da, b, &array);
120:   for (k=zs; k<zs+zm; k++) {
121:     for (j=ys; j<ys+ym; j++) {
122:       for (i=xs; i<xs+xm; i++) {
123:         for (d=0; d<dof; d++) {
124:           array[k][j][i][d] = 12 * PETSC_PI * PETSC_PI
125:                            * PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
126:                            * PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
127:                            * PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
128:                            * Hx * Hy * Hz;
129:         }
130:       }
131:     }
132:   }
133:   DMDAVecRestoreArrayDOFWrite(da, b, &array);
134:   VecAssemblyBegin(b);
135:   VecAssemblyEnd(b);

137:   /* force right hand side to be consistent for singular matrix */
138:   /* note this is really a hack, normally the model would provide you with a consistent right handside */

140:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
141:   MatNullSpaceRemove(nullspace,b);
142:   MatNullSpaceDestroy(&nullspace);
143:   return 0;
144: }

146: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac, void *ctx)
147: {
148:   PetscInt       dof,i,j,k,d,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
149:   PetscScalar    v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
150:   MatStencil     row, col[7];
151:   DM             da;
152:   MatNullSpace   nullspace;
153:   PetscBool      dump_mat = PETSC_FALSE, check_matis = PETSC_FALSE;

156:   KSPGetDM(ksp,&da);
157:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
158:   Hx      = 1.0 / (PetscReal)(mx);
159:   Hy      = 1.0 / (PetscReal)(my);
160:   Hz      = 1.0 / (PetscReal)(mz);
161:   HyHzdHx = Hy*Hz/Hx;
162:   HxHzdHy = Hx*Hz/Hy;
163:   HxHydHz = Hx*Hy/Hz;
164:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
165:   for (k=zs; k<zs+zm; k++) {
166:     for (j=ys; j<ys+ym; j++) {
167:       for (i=xs; i<xs+xm; i++) {
168:         for (d=0; d<dof; d++) {
169:           row.i = i; row.j = j; row.k = k; row.c = d;
170:           if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
171:             num = 0; numi=0; numj=0; numk=0;
172:             if (k!=0) {
173:               v[num]     = -HxHydHz;
174:               col[num].i = i;
175:               col[num].j = j;
176:               col[num].k = k-1;
177:               col[num].c = d;
178:               num++; numk++;
179:             }
180:             if (j!=0) {
181:               v[num]     = -HxHzdHy;
182:               col[num].i = i;
183:               col[num].j = j-1;
184:               col[num].k = k;
185:               col[num].c = d;
186:               num++; numj++;
187:               }
188:             if (i!=0) {
189:               v[num]     = -HyHzdHx;
190:               col[num].i = i-1;
191:               col[num].j = j;
192:               col[num].k = k;
193:               col[num].c = d;
194:               num++; numi++;
195:             }
196:             if (i!=mx-1) {
197:               v[num]     = -HyHzdHx;
198:               col[num].i = i+1;
199:               col[num].j = j;
200:               col[num].k = k;
201:               col[num].c = d;
202:               num++; numi++;
203:             }
204:             if (j!=my-1) {
205:               v[num]     = -HxHzdHy;
206:               col[num].i = i;
207:               col[num].j = j+1;
208:               col[num].k = k;
209:               col[num].c = d;
210:               num++; numj++;
211:             }
212:             if (k!=mz-1) {
213:               v[num]     = -HxHydHz;
214:               col[num].i = i;
215:               col[num].j = j;
216:               col[num].k = k+1;
217:               col[num].c = d;
218:               num++; numk++;
219:             }
220:             v[num]     = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
221:             col[num].i = i;   col[num].j = j;   col[num].k = k; col[num].c = d;
222:             num++;
223:             MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
224:           } else {
225:             v[0] = -HxHydHz;                          col[0].i = i;   col[0].j = j;   col[0].k = k-1; col[0].c = d;
226:             v[1] = -HxHzdHy;                          col[1].i = i;   col[1].j = j-1; col[1].k = k;   col[1].c = d;
227:             v[2] = -HyHzdHx;                          col[2].i = i-1; col[2].j = j;   col[2].k = k;   col[2].c = d;
228:             v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i;   col[3].j = j;   col[3].k = k;   col[3].c = d;
229:             v[4] = -HyHzdHx;                          col[4].i = i+1; col[4].j = j;   col[4].k = k;   col[4].c = d;
230:             v[5] = -HxHzdHy;                          col[5].i = i;   col[5].j = j+1; col[5].k = k;   col[5].c = d;
231:             v[6] = -HxHydHz;                          col[6].i = i;   col[6].j = j;   col[6].k = k+1; col[6].c = d;
232:             MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
233:           }
234:         }
235:       }
236:     }
237:   }
238:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
239:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
240:   PetscOptionsGetBool(NULL,NULL,"-dump_mat",&dump_mat,NULL);
241:   if (dump_mat) {
242:     Mat JJ;

244:     MatComputeOperator(jac,MATAIJ,&JJ);
245:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
246:     MatView(JJ,PETSC_VIEWER_STDOUT_WORLD);
247:     MatDestroy(&JJ);
248:   }
249:   MatViewFromOptions(jac,NULL,"-view_mat");
250:   PetscOptionsGetBool(NULL,NULL,"-check_matis",&check_matis,NULL);
251:   if (check_matis) {
252:     void      (*f)(void);
253:     Mat       J2;
254:     MatType   jtype;
255:     PetscReal nrm;

257:     MatGetType(jac,&jtype);
258:     MatConvert(jac,MATIS,MAT_INITIAL_MATRIX,&J2);
259:     MatViewFromOptions(J2,NULL,"-view_conv");
260:     MatConvert(J2,jtype,MAT_INPLACE_MATRIX,&J2);
261:     MatGetOperation(jac,MATOP_VIEW,&f);
262:     MatSetOperation(J2,MATOP_VIEW,f);
263:     MatSetDM(J2,da);
264:     MatViewFromOptions(J2,NULL,"-view_conv_assembled");
265:     MatAXPY(J2,-1.,jac,DIFFERENT_NONZERO_PATTERN);
266:     MatNorm(J2,NORM_FROBENIUS,&nrm);
267:     PetscPrintf(PETSC_COMM_WORLD,"Error MATIS %g\n",(double)nrm);
268:     MatViewFromOptions(J2,NULL,"-view_conv_err");
269:     MatDestroy(&J2);
270:   }
271:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
272:   MatSetNullSpace(J,nullspace);
273:   MatNullSpaceDestroy(&nullspace);
274:   return 0;
275: }

277: /*TEST

279:    build:
280:       requires: !complex !single

282:    test:
283:       args: -pc_type mg -pc_mg_type full -ksp_type fgmres -ksp_monitor_short -pc_mg_levels 3 -mg_coarse_pc_factor_shift_type nonzero -ksp_view

285:    test:
286:       suffix: 2
287:       nsize: 2
288:       args: -ksp_monitor_short -da_grid_x 50 -da_grid_y 50 -pc_type ksp -ksp_ksp_type cg -ksp_pc_type bjacobi -ksp_ksp_rtol 1e-1 -ksp_ksp_monitor -ksp_type pipefgmres -ksp_gmres_restart 5

290:    test:
291:       suffix: hyprestruct
292:       nsize: 3
293:       requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
294:       args: -ksp_type gmres -pc_type pfmg -dm_mat_type hyprestruct -ksp_monitor -da_refine 3

296: TEST*/