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;
37: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,d,dof;
38: PetscScalar Hx,Hy,Hz;
39: PetscScalar ****array;
40: Vec x,b,r;
41: Mat J;
43: PetscInitialize(&argc,&argv,(char*)0,help);
44: dof = 1;
45: PetscOptionsGetInt(NULL,NULL,"-da_dof",&dof,NULL);
46: KSPCreate(PETSC_COMM_WORLD,&ksp);
47: 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);
48: DMSetFromOptions(da);
49: DMSetUp(da);
50: DMDASetInterpolationType(da, DMDA_Q0);
52: KSPSetDM(ksp,da);
54: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
55: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
56: KSPSetFromOptions(ksp);
57: KSPSolve(ksp,NULL,NULL);
58: KSPGetSolution(ksp,&x);
59: KSPGetRhs(ksp,&b);
60: KSPGetOperators(ksp,NULL,&J);
61: VecDuplicate(b,&r);
63: MatMult(J,x,r);
64: VecAXPY(r,-1.0,b);
65: VecNorm(r,NORM_2,&norm);
66: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
68: DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
69: Hx = 1.0 / (PetscReal)(mx);
70: Hy = 1.0 / (PetscReal)(my);
71: Hz = 1.0 / (PetscReal)(mz);
72: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
73: DMDAVecGetArrayDOF(da, x, &array);
75: for (k=zs; k<zs+zm; k++) {
76: for (j=ys; j<ys+ym; j++) {
77: for (i=xs; i<xs+xm; i++) {
78: for (d=0; d<dof; d++) {
79: array[k][j][i][d] -=
80: PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
81: PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
82: PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
83: }
84: }
85: }
86: }
87: DMDAVecRestoreArrayDOF(da, x, &array);
88: VecAssemblyBegin(x);
89: VecAssemblyEnd(x);
91: VecNorm(x,NORM_INFINITY,&norm);
92: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)norm);
93: VecNorm(x,NORM_1,&norm);
94: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
95: VecNorm(x,NORM_2,&norm);
96: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
97: VecSum(x,&norm);
98: if (norm > 10000*PETSC_MACHINE_EPSILON) {
99: PetscPrintf(PETSC_COMM_WORLD,"Vector sum %g\n",(double)norm);
100: }
102: VecDestroy(&r);
103: KSPDestroy(&ksp);
104: DMDestroy(&da);
105: PetscFinalize();
106: return 0;
107: }
109: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
110: {
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: DMDAVecGetArrayDOFWrite(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: DMDAVecRestoreArrayDOFWrite(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: }
151: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac, void *ctx)
152: {
153: PetscInt dof,i,j,k,d,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
154: PetscScalar v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
155: MatStencil row, col[7];
156: DM da;
157: MatNullSpace nullspace;
158: PetscBool dump_mat = PETSC_FALSE, check_matis = PETSC_FALSE;
161: KSPGetDM(ksp,&da);
162: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,&dof,0,0,0,0,0);
163: Hx = 1.0 / (PetscReal)(mx);
164: Hy = 1.0 / (PetscReal)(my);
165: Hz = 1.0 / (PetscReal)(mz);
166: HyHzdHx = Hy*Hz/Hx;
167: HxHzdHy = Hx*Hz/Hy;
168: HxHydHz = Hx*Hy/Hz;
169: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
170: for (k=zs; k<zs+zm; k++) {
171: for (j=ys; j<ys+ym; j++) {
172: for (i=xs; i<xs+xm; i++) {
173: for (d=0; d<dof; d++) {
174: row.i = i; row.j = j; row.k = k; row.c = d;
175: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
176: num = 0; numi=0; numj=0; numk=0;
177: if (k!=0) {
178: v[num] = -HxHydHz;
179: col[num].i = i;
180: col[num].j = j;
181: col[num].k = k-1;
182: col[num].c = d;
183: num++; numk++;
184: }
185: if (j!=0) {
186: v[num] = -HxHzdHy;
187: col[num].i = i;
188: col[num].j = j-1;
189: col[num].k = k;
190: col[num].c = d;
191: num++; numj++;
192: }
193: if (i!=0) {
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++; numi++;
200: }
201: if (i!=mx-1) {
202: v[num] = -HyHzdHx;
203: col[num].i = i+1;
204: col[num].j = j;
205: col[num].k = k;
206: col[num].c = d;
207: num++; numi++;
208: }
209: if (j!=my-1) {
210: v[num] = -HxHzdHy;
211: col[num].i = i;
212: col[num].j = j+1;
213: col[num].k = k;
214: col[num].c = d;
215: num++; numj++;
216: }
217: if (k!=mz-1) {
218: v[num] = -HxHydHz;
219: col[num].i = i;
220: col[num].j = j;
221: col[num].k = k+1;
222: col[num].c = d;
223: num++; numk++;
224: }
225: v[num] = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
226: col[num].i = i; col[num].j = j; col[num].k = k; col[num].c = d;
227: num++;
228: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
229: } else {
230: v[0] = -HxHydHz; col[0].i = i; col[0].j = j; col[0].k = k-1; col[0].c = d;
231: v[1] = -HxHzdHy; col[1].i = i; col[1].j = j-1; col[1].k = k; col[1].c = d;
232: v[2] = -HyHzdHx; col[2].i = i-1; col[2].j = j; col[2].k = k; col[2].c = d;
233: v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i; col[3].j = j; col[3].k = k; col[3].c = d;
234: v[4] = -HyHzdHx; col[4].i = i+1; col[4].j = j; col[4].k = k; col[4].c = d;
235: v[5] = -HxHzdHy; col[5].i = i; col[5].j = j+1; col[5].k = k; col[5].c = d;
236: v[6] = -HxHydHz; col[6].i = i; col[6].j = j; col[6].k = k+1; col[6].c = d;
237: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
238: }
239: }
240: }
241: }
242: }
243: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
244: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
245: PetscOptionsGetBool(NULL,NULL,"-dump_mat",&dump_mat,NULL);
246: if (dump_mat) {
247: Mat JJ;
249: MatComputeOperator(jac,MATAIJ,&JJ);
250: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
251: MatView(JJ,PETSC_VIEWER_STDOUT_WORLD);
252: MatDestroy(&JJ);
253: }
254: MatViewFromOptions(jac,NULL,"-view_mat");
255: PetscOptionsGetBool(NULL,NULL,"-check_matis",&check_matis,NULL);
256: if (check_matis) {
257: void (*f)(void);
258: Mat J2;
259: MatType jtype;
260: PetscReal nrm;
262: MatGetType(jac,&jtype);
263: MatConvert(jac,MATIS,MAT_INITIAL_MATRIX,&J2);
264: MatViewFromOptions(J2,NULL,"-view_conv");
265: MatConvert(J2,jtype,MAT_INPLACE_MATRIX,&J2);
266: MatGetOperation(jac,MATOP_VIEW,&f);
267: MatSetOperation(J2,MATOP_VIEW,f);
268: MatSetDM(J2,da);
269: MatViewFromOptions(J2,NULL,"-view_conv_assembled");
270: MatAXPY(J2,-1.,jac,DIFFERENT_NONZERO_PATTERN);
271: MatNorm(J2,NORM_FROBENIUS,&nrm);
272: PetscPrintf(PETSC_COMM_WORLD,"Error MATIS %g\n",(double)nrm);
273: MatViewFromOptions(J2,NULL,"-view_conv_err");
274: MatDestroy(&J2);
275: }
276: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
277: MatSetNullSpace(J,nullspace);
278: MatNullSpaceDestroy(&nullspace);
279: return 0;
280: }
282: /*TEST
284: build:
285: requires: !complex !single
287: test:
288: 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
290: test:
291: suffix: 2
292: nsize: 2
293: 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
295: test:
296: suffix: hyprestruct
297: nsize: 3
298: requires: hypre !defined(PETSC_HAVE_HYPRE_DEVICE)
299: args: -ksp_type gmres -pc_type pfmg -dm_mat_type hyprestruct -ksp_monitor -da_refine 3
301: TEST*/