Actual source code: ex34.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
```  1: /*T
2:    Concepts: KSP^solving a system of linear equations
3:    Concepts: KSP^Laplacian, 3d
4:    Processors: n
5: T*/

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

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

14: with pure Neumann boundary conditions

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

18: The functions are cell-centered

20: This uses multigrid to solve the linear system

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

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

27: #include <petscdm.h>
28: #include <petscdmda.h>
29: #include <petscksp.h>

31: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
32: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);

34: int main(int argc,char **argv)
35: {
36:   KSP            ksp;
37:   DM             da;
38:   PetscReal      norm;
40:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,d,dof;
41:   PetscScalar    Hx,Hy,Hz;
42:   PetscScalar    ****array;
43:   Vec            x,b,r;
44:   Mat            J;

46:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
47:   dof  = 1;
48:   PetscOptionsGetInt(NULL,NULL,"-da_dof",&dof,NULL);
49:   KSPCreate(PETSC_COMM_WORLD,&ksp);
50:   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);
51:   DMSetFromOptions(da);
52:   DMSetUp(da);
53:   DMDASetInterpolationType(da, DMDA_Q0);

55:   KSPSetDM(ksp,da);

57:   KSPSetComputeRHS(ksp,ComputeRHS,NULL);
58:   KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
59:   KSPSetFromOptions(ksp);
60:   KSPSolve(ksp,NULL,NULL);
61:   KSPGetSolution(ksp,&x);
62:   KSPGetRhs(ksp,&b);
63:   KSPGetOperators(ksp,NULL,&J);
64:   VecDuplicate(b,&r);

66:   MatMult(J,x,r);
67:   VecAXPY(r,-1.0,b);
68:   VecNorm(r,NORM_2,&norm);
69:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

71:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
72:   Hx   = 1.0 / (PetscReal)(mx);
73:   Hy   = 1.0 / (PetscReal)(my);
74:   Hz   = 1.0 / (PetscReal)(mz);
75:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
76:   DMDAVecGetArrayDOF(da, x, &array);

78:   for (k=zs; k<zs+zm; k++) {
79:     for (j=ys; j<ys+ym; j++) {
80:       for (i=xs; i<xs+xm; i++) {
81:         for (d=0; d<dof; d++) {
82:           array[k][j][i][d] -=
83:             PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
84:             PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
85:             PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
86:         }
87:       }
88:     }
89:   }
90:   DMDAVecRestoreArrayDOF(da, x, &array);
91:   VecAssemblyBegin(x);
92:   VecAssemblyEnd(x);

94:   VecNorm(x,NORM_INFINITY,&norm);
95:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)norm);
96:   VecNorm(x,NORM_1,&norm);
97:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
98:   VecNorm(x,NORM_2,&norm);
99:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));

101:   VecDestroy(&r);
102:   KSPDestroy(&ksp);
103:   DMDestroy(&da);
104:   PetscFinalize();
105:   return ierr;
106: }

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

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

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

145:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
146:   MatNullSpaceRemove(nullspace,b);
147:   MatNullSpaceDestroy(&nullspace);
148:   return(0);
149: }

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

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

251:     MatComputeOperator(jac,MATAIJ,&JJ);
252:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
253:     MatView(JJ,PETSC_VIEWER_STDOUT_WORLD);
254:     MatDestroy(&JJ);
255:   }
256:   MatViewFromOptions(jac,NULL,"-view_mat");
257:   PetscOptionsGetBool(NULL,NULL,"-check_matis",&check_matis,NULL);
258:   if (check_matis) {
259:     void      (*f)(void);
260:     Mat       J2;
261:     MatType   jtype;
262:     PetscReal nrm;

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

286: /*TEST

288:    build:
289:       requires: !complex !single

291:    test:
292:       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

294:    test:
295:       suffix: 2
296:       nsize: 2
297:       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

299:    test:
300:       suffix: hyprestruct
301:       nsize: 3
302:       requires: hypre
303:       args: -ksp_type gmres -pc_type pfmg -dm_mat_type hyprestruct -ksp_monitor -da_refine 3

305: TEST*/
```