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*/