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: }

```