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;
39: PetscInitialize(&argc, &argv, (char *)0, help);
40: dof = 1;
41: PetscOptionsGetInt(NULL, NULL, "-da_dof", &dof, NULL);
42: KSPCreate(PETSC_COMM_WORLD, &ksp);
43: 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);
44: DMSetFromOptions(da);
45: DMSetUp(da);
46: DMDASetInterpolationType(da, DMDA_Q0);
48: KSPSetDM(ksp, da);
50: KSPSetComputeRHS(ksp, ComputeRHS, NULL);
51: KSPSetComputeOperators(ksp, ComputeMatrix, NULL);
52: KSPSetFromOptions(ksp);
53: KSPSolve(ksp, NULL, NULL);
54: KSPGetSolution(ksp, &x);
55: KSPGetRhs(ksp, &b);
56: KSPGetOperators(ksp, NULL, &J);
57: VecDuplicate(b, &r);
59: MatMult(J, x, r);
60: VecAXPY(r, -1.0, b);
61: VecNorm(r, NORM_2, &norm);
62: PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm);
64: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0);
65: Hx = 1.0 / (PetscReal)(mx);
66: Hy = 1.0 / (PetscReal)(my);
67: Hz = 1.0 / (PetscReal)(mz);
68: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
69: DMDAVecGetArrayDOF(da, x, &array);
71: for (k = zs; k < zs + zm; k++) {
72: for (j = ys; j < ys + ym; j++) {
73: for (i = xs; i < xs + xm; i++) {
74: for (d = 0; d < dof; d++) array[k][j][i][d] -= PetscCosScalar(2 * PETSC_PI * (((PetscReal)i + 0.5) * Hx)) * PetscCosScalar(2 * PETSC_PI * (((PetscReal)j + 0.5) * Hy)) * PetscCosScalar(2 * PETSC_PI * (((PetscReal)k + 0.5) * Hz));
75: }
76: }
77: }
78: DMDAVecRestoreArrayDOF(da, x, &array);
79: VecAssemblyBegin(x);
80: VecAssemblyEnd(x);
82: VecNorm(x, NORM_INFINITY, &norm);
83: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)norm);
84: VecNorm(x, NORM_1, &norm);
85: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)(norm / ((PetscReal)(mx) * (PetscReal)(my) * (PetscReal)(mz))));
86: VecNorm(x, NORM_2, &norm);
87: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)(norm / ((PetscReal)(mx) * (PetscReal)(my) * (PetscReal)(mz))));
89: VecDestroy(&r);
90: KSPDestroy(&ksp);
91: DMDestroy(&da);
92: PetscFinalize();
93: return 0;
94: }
96: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
97: {
98: PetscInt d, dof, i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs;
99: PetscScalar Hx, Hy, Hz;
100: PetscScalar ****array;
101: DM da;
102: MatNullSpace nullspace;
105: KSPGetDM(ksp, &da);
106: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
107: Hx = 1.0 / (PetscReal)(mx);
108: Hy = 1.0 / (PetscReal)(my);
109: Hz = 1.0 / (PetscReal)(mz);
110: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
111: DMDAVecGetArrayDOFWrite(da, b, &array);
112: for (k = zs; k < zs + zm; k++) {
113: for (j = ys; j < ys + ym; j++) {
114: for (i = xs; i < xs + xm; i++) {
115: for (d = 0; d < dof; d++) {
116: array[k][j][i][d] = 12 * PETSC_PI * PETSC_PI * PetscCosScalar(2 * PETSC_PI * (((PetscReal)i + 0.5) * Hx)) * PetscCosScalar(2 * PETSC_PI * (((PetscReal)j + 0.5) * Hy)) * PetscCosScalar(2 * PETSC_PI * (((PetscReal)k + 0.5) * Hz)) * Hx * Hy * Hz;
117: }
118: }
119: }
120: }
121: DMDAVecRestoreArrayDOFWrite(da, b, &array);
122: VecAssemblyBegin(b);
123: VecAssemblyEnd(b);
125: /* force right hand side to be consistent for singular matrix */
126: /* note this is really a hack, normally the model would provide you with a consistent right handside */
128: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
129: MatNullSpaceRemove(nullspace, b);
130: MatNullSpaceDestroy(&nullspace);
131: return 0;
132: }
134: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
135: {
136: PetscInt dof, i, j, k, d, mx, my, mz, xm, ym, zm, xs, ys, zs, num, numi, numj, numk;
137: PetscScalar v[7], Hx, Hy, Hz, HyHzdHx, HxHzdHy, HxHydHz;
138: MatStencil row, col[7];
139: DM da;
140: MatNullSpace nullspace;
141: PetscBool dump_mat = PETSC_FALSE, check_matis = PETSC_FALSE;
144: KSPGetDM(ksp, &da);
145: DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
146: Hx = 1.0 / (PetscReal)(mx);
147: Hy = 1.0 / (PetscReal)(my);
148: Hz = 1.0 / (PetscReal)(mz);
149: HyHzdHx = Hy * Hz / Hx;
150: HxHzdHy = Hx * Hz / Hy;
151: HxHydHz = Hx * Hy / Hz;
152: DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
153: for (k = zs; k < zs + zm; k++) {
154: for (j = ys; j < ys + ym; j++) {
155: for (i = xs; i < xs + xm; i++) {
156: for (d = 0; d < dof; d++) {
157: row.i = i;
158: row.j = j;
159: row.k = k;
160: row.c = d;
161: if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || k == mz - 1) {
162: num = 0;
163: numi = 0;
164: numj = 0;
165: numk = 0;
166: if (k != 0) {
167: v[num] = -HxHydHz;
168: col[num].i = i;
169: col[num].j = j;
170: col[num].k = k - 1;
171: col[num].c = d;
172: num++;
173: numk++;
174: }
175: if (j != 0) {
176: v[num] = -HxHzdHy;
177: col[num].i = i;
178: col[num].j = j - 1;
179: col[num].k = k;
180: col[num].c = d;
181: num++;
182: numj++;
183: }
184: if (i != 0) {
185: v[num] = -HyHzdHx;
186: col[num].i = i - 1;
187: col[num].j = j;
188: col[num].k = k;
189: col[num].c = d;
190: num++;
191: numi++;
192: }
193: if (i != mx - 1) {
194: v[num] = -HyHzdHx;
195: col[num].i = i + 1;
196: col[num].j = j;
197: col[num].k = k;
198: col[num].c = d;
199: num++;
200: numi++;
201: }
202: if (j != my - 1) {
203: v[num] = -HxHzdHy;
204: col[num].i = i;
205: col[num].j = j + 1;
206: col[num].k = k;
207: col[num].c = d;
208: num++;
209: numj++;
210: }
211: if (k != mz - 1) {
212: v[num] = -HxHydHz;
213: col[num].i = i;
214: col[num].j = j;
215: col[num].k = k + 1;
216: col[num].c = d;
217: num++;
218: numk++;
219: }
220: v[num] = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
221: col[num].i = i;
222: col[num].j = j;
223: col[num].k = k;
224: col[num].c = d;
225: num++;
226: MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES);
227: } else {
228: v[0] = -HxHydHz;
229: col[0].i = i;
230: col[0].j = j;
231: col[0].k = k - 1;
232: col[0].c = d;
233: v[1] = -HxHzdHy;
234: col[1].i = i;
235: col[1].j = j - 1;
236: col[1].k = k;
237: col[1].c = d;
238: v[2] = -HyHzdHx;
239: col[2].i = i - 1;
240: col[2].j = j;
241: col[2].k = k;
242: col[2].c = d;
243: v[3] = 2.0 * (HyHzdHx + HxHzdHy + HxHydHz);
244: col[3].i = i;
245: col[3].j = j;
246: col[3].k = k;
247: col[3].c = d;
248: v[4] = -HyHzdHx;
249: col[4].i = i + 1;
250: col[4].j = j;
251: col[4].k = k;
252: col[4].c = d;
253: v[5] = -HxHzdHy;
254: col[5].i = i;
255: col[5].j = j + 1;
256: col[5].k = k;
257: col[5].c = d;
258: v[6] = -HxHydHz;
259: col[6].i = i;
260: col[6].j = j;
261: col[6].k = k + 1;
262: col[6].c = d;
263: MatSetValuesStencil(jac, 1, &row, 7, col, v, INSERT_VALUES);
264: }
265: }
266: }
267: }
268: }
269: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
270: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
271: PetscOptionsGetBool(NULL, NULL, "-dump_mat", &dump_mat, NULL);
272: if (dump_mat) {
273: Mat JJ;
275: MatComputeOperator(jac, MATAIJ, &JJ);
276: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_MATLAB);
277: MatView(JJ, PETSC_VIEWER_STDOUT_WORLD);
278: MatDestroy(&JJ);
279: }
280: MatViewFromOptions(jac, NULL, "-view_mat");
281: PetscOptionsGetBool(NULL, NULL, "-check_matis", &check_matis, NULL);
282: if (check_matis) {
283: void (*f)(void);
284: Mat J2;
285: MatType jtype;
286: PetscReal nrm;
288: MatGetType(jac, &jtype);
289: MatConvert(jac, MATIS, MAT_INITIAL_MATRIX, &J2);
290: MatViewFromOptions(J2, NULL, "-view_conv");
291: MatConvert(J2, jtype, MAT_INPLACE_MATRIX, &J2);
292: MatGetOperation(jac, MATOP_VIEW, &f);
293: MatSetOperation(J2, MATOP_VIEW, f);
294: MatSetDM(J2, da);
295: MatViewFromOptions(J2, NULL, "-view_conv_assembled");
296: MatAXPY(J2, -1., jac, DIFFERENT_NONZERO_PATTERN);
297: MatNorm(J2, NORM_FROBENIUS, &nrm);
298: PetscPrintf(PETSC_COMM_WORLD, "Error MATIS %g\n", (double)nrm);
299: MatViewFromOptions(J2, NULL, "-view_conv_err");
300: MatDestroy(&J2);
301: }
302: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
303: MatSetNullSpace(J, nullspace);
304: MatNullSpaceDestroy(&nullspace);
305: return 0;
306: }
308: /*TEST
310: build:
311: requires: !complex !single
313: test:
314: 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
316: test:
317: suffix: 2
318: nsize: 2
319: 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
321: test:
322: suffix: hyprestruct
323: nsize: 3
324: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
325: args: -ksp_type gmres -pc_type pfmg -dm_mat_type hyprestruct -ksp_monitor -da_refine 3
327: TEST*/