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