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 petscda.h
26: #include petscksp.h
27: #include petscmg.h
28: #include petscdmmg.h
33: typedef enum {DIRICHLET, NEUMANN} BCType;
35: typedef struct {
36: BCType bcType;
37: } UserContext;
41: int main(int argc,char **argv)
42: {
43: DMMG *dmmg;
44: DA da;
45: UserContext user;
46: PetscReal norm;
47: const char *bcTypes[2] = {"dirichlet","neumann"};
49: PetscInt l,bc;
51: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
52: PetscScalar Hx,Hy,Hz;
53: PetscScalar ***array;
56: PetscInitialize(&argc,&argv,(char *)0,help);
57:
58: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
59: DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
60: DASetInterpolationType(da, DA_Q0);
62: DMMGSetDM(dmmg,(DM)da);
63: DADestroy(da);
64: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
65: DMMGSetUser(dmmg,l,&user);
66: }
67:
68: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DMMG");
69: bc = (PetscInt)NEUMANN;
70: PetscOptionsEList("-bc_type","Type of boundary condition","ex34.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
71: user.bcType = (BCType)bc;
72: PetscOptionsEnd();
73:
74: DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
75: if (user.bcType == NEUMANN) {
76: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
77: }
79: DMMGSolve(dmmg);
80:
81: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
82: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
83: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
84: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
85:
86: DAGetInfo(DMMGGetDA(dmmg), 0, &mx, &my, &mz, 0,0,0,0,0,0,0);
87: Hx = 1.0 / (PetscReal)(mx);
88: Hy = 1.0 / (PetscReal)(my);
89: Hz = 1.0 / (PetscReal)(mz);
90: DAGetCorners(DMMGGetDA(dmmg),&xs,&ys,&zs,&xm,&ym,&zm);
91: DAVecGetArray(DMMGGetDA(dmmg), DMMGGetx(dmmg), &array);
93: for (k=zs; k<zs+zm; k++){
94: for (j=ys; j<ys+ym; j++){
95: for(i=xs; i<xs+xm; i++){
96: array[k][j][i] -=
97: PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
98: PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
99: PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
100: }
101: }
102: }
103: DAVecRestoreArray(DMMGGetDA(dmmg), DMMGGetx(dmmg), &array);
104: VecAssemblyBegin(DMMGGetx(dmmg));
105: VecAssemblyEnd(DMMGGetx(dmmg));
107: VecNorm(DMMGGetx(dmmg),NORM_INFINITY,&norm);
108: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm);
109: VecNorm(DMMGGetx(dmmg),NORM_1,&norm);
110: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
111: VecNorm(DMMGGetx(dmmg),NORM_2,&norm);
112: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
114: DMMGDestroy(dmmg);
115: PetscFinalize();
116: return 0;
117: }
121: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
122: {
123: DA da = (DA)dmmg->dm;
124: UserContext *user = (UserContext *) dmmg->user;
126: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
127: PetscScalar Hx,Hy,Hz;
128: PetscScalar ***array;
131: DAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0);
132: Hx = 1.0 / (PetscReal)(mx);
133: Hy = 1.0 / (PetscReal)(my);
134: Hz = 1.0 / (PetscReal)(mz);
135: DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
136: DAVecGetArray(da, b, &array);
137: for (k=zs; k<zs+zm; k++){
138: for (j=ys; j<ys+ym; j++){
139: for(i=xs; i<xs+xm; i++){
140: array[k][j][i] = 12*PETSC_PI*PETSC_PI
141: *PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
142: *PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
143: *PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
144: *Hx*Hy*Hz;
145: }
146: }
147: }
148: DAVecRestoreArray(da, b, &array);
149: VecAssemblyBegin(b);
150: VecAssemblyEnd(b);
152: /* force right hand side to be consistent for singular matrix */
153: /* note this is really a hack, normally the model would provide you with a consistent right handside */
154: if (user->bcType == NEUMANN)
155: {
156: MatNullSpace nullspace;
157:
158: KSPGetNullSpace(dmmg->ksp,&nullspace);
159: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
160: }
162: return(0);
163: }
165:
168: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J,Mat jac)
169: {
170: DA da = (DA) dmmg->dm;
171: UserContext *user = (UserContext *) dmmg->user;
173: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
174: PetscScalar v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
175: MatStencil row, col[7];
178: DAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0);
179: Hx = 1.0 / (PetscReal)(mx);
180: Hy = 1.0 / (PetscReal)(my);
181: Hz = 1.0 / (PetscReal)(mz);
182: HyHzdHx = Hy*Hz/Hx;
183: HxHzdHy = Hx*Hz/Hy;
184: HxHydHz = Hx*Hy/Hz;
185: DAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
186: for (k=zs; k<zs+zm; k++)
187: {
188: for (j=ys; j<ys+ym; j++)
189: {
190: for(i=xs; i<xs+xm; i++)
191: {
192: row.i = i; row.j = j; row.k = k;
193: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1)
194: {
195: if (user->bcType == DIRICHLET)
196: {
197: SETERRQ(PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
198: v[0] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz);
199: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
200: }
201: else if (user->bcType == NEUMANN)
202: {
203: num = 0; numi=0; numj=0; numk=0;
204: if (k!=0)
205: {
206: v[num] = -HxHydHz;
207: col[num].i = i;
208: col[num].j = j;
209: col[num].k = k-1;
210: num++; numk++;
211: }
212: if (j!=0)
213: {
214: v[num] = -HxHzdHy;
215: col[num].i = i;
216: col[num].j = j-1;
217: col[num].k = k;
218: num++; numj++;
219: }
220: if (i!=0)
221: {
222: v[num] = -HyHzdHx;
223: col[num].i = i-1;
224: col[num].j = j;
225: col[num].k = k;
226: num++; numi++;
227: }
228: if (i!=mx-1)
229: {
230: v[num] = -HyHzdHx;
231: col[num].i = i+1;
232: col[num].j = j;
233: col[num].k = k;
234: num++; numi++;
235: }
236: if (j!=my-1)
237: {
238: v[num] = -HxHzdHy;
239: col[num].i = i;
240: col[num].j = j+1;
241: col[num].k = k;
242: num++; numj++;
243: }
244: if (k!=mz-1)
245: {
246: v[num] = -HxHydHz;
247: col[num].i = i;
248: col[num].j = j;
249: col[num].k = k+1;
250: num++; numk++;
251: }
252: v[num] = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
253: col[num].i = i; col[num].j = j; col[num].k = k;
254: num++;
255: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
256: }
257: }
258: else
259: {
260: v[0] = -HxHydHz; col[0].i = i; col[0].j = j; col[0].k = k-1;
261: v[1] = -HxHzdHy; col[1].i = i; col[1].j = j-1; col[1].k = k;
262: v[2] = -HyHzdHx; col[2].i = i-1; col[2].j = j; col[2].k = k;
263: v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i; col[3].j = j; col[3].k = k;
264: v[4] = -HyHzdHx; col[4].i = i+1; col[4].j = j; col[4].k = k;
265: v[5] = -HxHzdHy; col[5].i = i; col[5].j = j+1; col[5].k = k;
266: v[6] = -HxHydHz; col[6].i = i; col[6].j = j; col[6].k = k+1;
267: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
268: }
269: }
270: }
271: }
272: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
273: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
275: return(0);
276: }