Actual source code: ex34.c
petsc-3.3-p7 2013-05-11
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 <petscdmda.h>
26: #include <petscksp.h>
27: #include <petscpcmg.h>
29: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
30: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
32: typedef enum {DIRICHLET, NEUMANN} BCType;
34: typedef struct {
35: BCType bcType;
36: } UserContext;
40: int main(int argc,char **argv)
41: {
42: KSP ksp;
43: DM da;
44: UserContext user;
45: PetscReal norm;
46: const char *bcTypes[2] = {"dirichlet","neumann"};
48: PetscInt bc;
50: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
51: PetscScalar Hx,Hy,Hz;
52: PetscScalar ***array;
53: Vec x,b,r;
54: Mat J;
56: PetscInitialize(&argc,&argv,(char *)0,help);
57:
58: KSPCreate(PETSC_COMM_WORLD,&ksp);
59: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,12,12,12,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
60: DMDASetInterpolationType(da, DMDA_Q0);
62: KSPSetDM(ksp,da);
63:
64: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
65: bc = (PetscInt)NEUMANN;
66: PetscOptionsEList("-bc_type","Type of boundary condition","ex34.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
67: user.bcType = (BCType)bc;
68: PetscOptionsEnd();
69:
70: KSPSetComputeRHS(ksp,ComputeRHS,&user);
71: KSPSetComputeOperators(ksp,ComputeMatrix,&user);
72: KSPSetFromOptions(ksp);
73: KSPSolve(ksp,PETSC_NULL,PETSC_NULL);
74: KSPGetSolution(ksp,&x);
75: KSPGetRhs(ksp,&b);
76: KSPGetOperators(ksp,PETSC_NULL,&J,PETSC_NULL);
77: VecDuplicate(b,&r);
79: MatMult(J,x,r);
80: VecAXPY(r,-1.0,b);
81: VecNorm(r,NORM_2,&norm);
82: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
83:
84: DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
85: Hx = 1.0 / (PetscReal)(mx);
86: Hy = 1.0 / (PetscReal)(my);
87: Hz = 1.0 / (PetscReal)(mz);
88: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
89: DMDAVecGetArray(da, x, &array);
91: for (k=zs; k<zs+zm; k++){
92: for (j=ys; j<ys+ym; j++){
93: for(i=xs; i<xs+xm; i++){
94: array[k][j][i] -=
95: PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
96: PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
97: PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
98: }
99: }
100: }
101: DMDAVecRestoreArray(da, x, &array);
102: VecAssemblyBegin(x);
103: VecAssemblyEnd(x);
105: VecNorm(x,NORM_INFINITY,&norm);
106: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm);
107: VecNorm(x,NORM_1,&norm);
108: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
109: VecNorm(x,NORM_2,&norm);
110: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
112: VecDestroy(&r);
113: KSPDestroy(&ksp);
114: DMDestroy(&da);
115: PetscFinalize();
116: return 0;
117: }
121: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
122: {
123: UserContext *user = (UserContext*)ctx;
125: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
126: PetscScalar Hx,Hy,Hz;
127: PetscScalar ***array;
128: DM da;
131: KSPGetDM(ksp,&da);
132: DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
133: Hx = 1.0 / (PetscReal)(mx);
134: Hy = 1.0 / (PetscReal)(my);
135: Hz = 1.0 / (PetscReal)(mz);
136: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
137: DMDAVecGetArray(da, b, &array);
138: for (k=zs; k<zs+zm; k++){
139: for (j=ys; j<ys+ym; j++){
140: for(i=xs; i<xs+xm; i++){
141: array[k][j][i] = 12*PETSC_PI*PETSC_PI
142: *PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
143: *PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
144: *PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
145: *Hx*Hy*Hz;
146: }
147: }
148: }
149: DMDAVecRestoreArray(da, b, &array);
150: VecAssemblyBegin(b);
151: VecAssemblyEnd(b);
153: /* force right hand side to be consistent for singular matrix */
154: /* note this is really a hack, normally the model would provide you with a consistent right handside */
155: if (user->bcType == NEUMANN) {
156: MatNullSpace nullspace;
158: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
159: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
160: MatNullSpaceDestroy(&nullspace);
161: }
162: return(0);
163: }
165:
168: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac,MatStructure *str, void *ctx)
169: {
170: UserContext *user = (UserContext*)ctx;
172: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
173: PetscScalar v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
174: MatStencil row, col[7];
175: DM da;
178: KSPGetDM(ksp,&da);
179: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
180: Hx = 1.0 / (PetscReal)(mx);
181: Hy = 1.0 / (PetscReal)(my);
182: Hz = 1.0 / (PetscReal)(mz);
183: HyHzdHx = Hy*Hz/Hx;
184: HxHzdHy = Hx*Hz/Hy;
185: HxHydHz = Hx*Hy/Hz;
186: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
187: for (k=zs; k<zs+zm; k++)
188: {
189: for (j=ys; j<ys+ym; j++)
190: {
191: for(i=xs; i<xs+xm; i++)
192: {
193: row.i = i; row.j = j; row.k = k;
194: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1)
195: {
196: if (user->bcType == DIRICHLET)
197: {
198: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
199: v[0] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz);
200: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
201: }
202: else if (user->bcType == NEUMANN)
203: {
204: num = 0; numi=0; numj=0; numk=0;
205: if (k!=0)
206: {
207: v[num] = -HxHydHz;
208: col[num].i = i;
209: col[num].j = j;
210: col[num].k = k-1;
211: num++; numk++;
212: }
213: if (j!=0)
214: {
215: v[num] = -HxHzdHy;
216: col[num].i = i;
217: col[num].j = j-1;
218: col[num].k = k;
219: num++; numj++;
220: }
221: if (i!=0)
222: {
223: v[num] = -HyHzdHx;
224: col[num].i = i-1;
225: col[num].j = j;
226: col[num].k = k;
227: num++; numi++;
228: }
229: if (i!=mx-1)
230: {
231: v[num] = -HyHzdHx;
232: col[num].i = i+1;
233: col[num].j = j;
234: col[num].k = k;
235: num++; numi++;
236: }
237: if (j!=my-1)
238: {
239: v[num] = -HxHzdHy;
240: col[num].i = i;
241: col[num].j = j+1;
242: col[num].k = k;
243: num++; numj++;
244: }
245: if (k!=mz-1)
246: {
247: v[num] = -HxHydHz;
248: col[num].i = i;
249: col[num].j = j;
250: col[num].k = k+1;
251: num++; numk++;
252: }
253: v[num] = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
254: col[num].i = i; col[num].j = j; col[num].k = k;
255: num++;
256: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
257: }
258: }
259: else
260: {
261: v[0] = -HxHydHz; col[0].i = i; col[0].j = j; col[0].k = k-1;
262: v[1] = -HxHzdHy; col[1].i = i; col[1].j = j-1; col[1].k = k;
263: v[2] = -HyHzdHx; col[2].i = i-1; col[2].j = j; col[2].k = k;
264: v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i; col[3].j = j; col[3].k = k;
265: v[4] = -HyHzdHx; col[4].i = i+1; col[4].j = j; col[4].k = k;
266: v[5] = -HxHzdHy; col[5].i = i; col[5].j = j+1; col[5].k = k;
267: v[6] = -HxHydHz; col[6].i = i; col[6].j = j; col[6].k = k+1;
268: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
269: }
270: }
271: }
272: }
273: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
274: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
275: if (user->bcType == NEUMANN) {
276: MatNullSpace nullspace;
278: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
279: MatSetNullSpace(jac,nullspace);
280: MatNullSpaceDestroy(&nullspace);
281: }
282: return(0);
283: }