Actual source code: ex34.c

petsc-3.4.4 2014-03-13
  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>

 28: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
 29: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);

 31: typedef enum {DIRICHLET, NEUMANN} BCType;

 33: typedef struct {
 34:   BCType bcType;
 35: } UserContext;

 39: int main(int argc,char **argv)
 40: {
 41:   KSP            ksp;
 42:   DM             da;
 43:   UserContext    user;
 44:   PetscReal      norm;
 45:   const char     *bcTypes[2] = {"dirichlet","neumann"};
 47:   PetscInt       bc;

 49:   PetscInt    i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
 50:   PetscScalar Hx,Hy,Hz;
 51:   PetscScalar ***array;
 52:   Vec         x,b,r;
 53:   Mat         J;

 55:   PetscInitialize(&argc,&argv,(char*)0,help);

 57:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 58:   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);
 59:   DMDASetInterpolationType(da, DMDA_Q0);

 61:   KSPSetDM(ksp,da);

 63:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
 64:   bc          = (PetscInt)NEUMANN;
 65:   PetscOptionsEList("-bc_type","Type of boundary condition","ex34.c",bcTypes,2,bcTypes[0],&bc,NULL);
 66:   user.bcType = (BCType)bc;
 67:   PetscOptionsEnd();

 69:   KSPSetComputeRHS(ksp,ComputeRHS,&user);
 70:   KSPSetComputeOperators(ksp,ComputeMatrix,&user);
 71:   KSPSetFromOptions(ksp);
 72:   KSPSolve(ksp,NULL,NULL);
 73:   KSPGetSolution(ksp,&x);
 74:   KSPGetRhs(ksp,&b);
 75:   KSPGetOperators(ksp,NULL,&J,NULL);
 76:   VecDuplicate(b,&r);

 78:   MatMult(J,x,r);
 79:   VecAXPY(r,-1.0,b);
 80:   VecNorm(r,NORM_2,&norm);
 81:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);

 83:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
 84:   Hx   = 1.0 / (PetscReal)(mx);
 85:   Hy   = 1.0 / (PetscReal)(my);
 86:   Hz   = 1.0 / (PetscReal)(mz);
 87:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
 88:   DMDAVecGetArray(da, x, &array);

 90:   for (k=zs; k<zs+zm; k++) {
 91:     for (j=ys; j<ys+ym; j++) {
 92:       for (i=xs; i<xs+xm; i++) {
 93:         array[k][j][i] -=
 94:           PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
 95:           PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
 96:           PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
 97:       }
 98:     }
 99:   }
100:   DMDAVecRestoreArray(da, x, &array);
101:   VecAssemblyBegin(x);
102:   VecAssemblyEnd(x);

104:   VecNorm(x,NORM_INFINITY,&norm);
105:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm);
106:   VecNorm(x,NORM_1,&norm);
107:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
108:   VecNorm(x,NORM_2,&norm);
109:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));

111:   VecDestroy(&r);
112:   KSPDestroy(&ksp);
113:   DMDestroy(&da);
114:   PetscFinalize();
115:   return 0;
116: }

120: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
121: {
122:   UserContext    *user = (UserContext*)ctx;
124:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
125:   PetscScalar    Hx,Hy,Hz;
126:   PetscScalar    ***array;
127:   DM             da;

130:   KSPGetDM(ksp,&da);
131:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,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:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
136:   DMDAVecGetArray(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:   DMDAVecRestoreArray(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:     MatNullSpace nullspace;

157:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
158:     MatNullSpaceRemove(nullspace,b,NULL);
159:     MatNullSpaceDestroy(&nullspace);
160:   }
161:   return(0);
162: }


167: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac,MatStructure *str, void *ctx)
168: {
169:   UserContext    *user = (UserContext*)ctx;
171:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
172:   PetscScalar    v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
173:   MatStencil     row, col[7];
174:   DM             da;

177:   KSPGetDM(ksp,&da);
178:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,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:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
186:   for (k=zs; k<zs+zm; k++) {
187:     for (j=ys; j<ys+ym; j++) {
188:       for (i=xs; i<xs+xm; i++) {
189:         row.i = i; row.j = j; row.k = k;
190:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
191:           if (user->bcType == DIRICHLET) {
192:             SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
193:             v[0] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz);
194:             MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
195:           } else if (user->bcType == NEUMANN) {
196:             num = 0; numi=0; numj=0; numk=0;
197:             if (k!=0) {
198:               v[num]     = -HxHydHz;
199:               col[num].i = i;
200:               col[num].j = j;
201:               col[num].k = k-1;
202:               num++; numk++;
203:             }
204:             if (j!=0) {
205:               v[num]     = -HxHzdHy;
206:               col[num].i = i;
207:               col[num].j = j-1;
208:               col[num].k = k;
209:               num++; numj++;
210:             }
211:             if (i!=0) {
212:               v[num]     = -HyHzdHx;
213:               col[num].i = i-1;
214:               col[num].j = j;
215:               col[num].k = k;
216:               num++; numi++;
217:             }
218:             if (i!=mx-1) {
219:               v[num]     = -HyHzdHx;
220:               col[num].i = i+1;
221:               col[num].j = j;
222:               col[num].k = k;
223:               num++; numi++;
224:             }
225:             if (j!=my-1) {
226:               v[num]     = -HxHzdHy;
227:               col[num].i = i;
228:               col[num].j = j+1;
229:               col[num].k = k;
230:               num++; numj++;
231:             }
232:             if (k!=mz-1) {
233:               v[num]     = -HxHydHz;
234:               col[num].i = i;
235:               col[num].j = j;
236:               col[num].k = k+1;
237:               num++; numk++;
238:             }
239:             v[num]     = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
240:             col[num].i = i;   col[num].j = j;   col[num].k = k;
241:             num++;
242:             MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
243:           }
244:         } else {
245:           v[0] = -HxHydHz;                          col[0].i = i;   col[0].j = j;   col[0].k = k-1;
246:           v[1] = -HxHzdHy;                          col[1].i = i;   col[1].j = j-1; col[1].k = k;
247:           v[2] = -HyHzdHx;                          col[2].i = i-1; col[2].j = j;   col[2].k = k;
248:           v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i;   col[3].j = j;   col[3].k = k;
249:           v[4] = -HyHzdHx;                          col[4].i = i+1; col[4].j = j;   col[4].k = k;
250:           v[5] = -HxHzdHy;                          col[5].i = i;   col[5].j = j+1; col[5].k = k;
251:           v[6] = -HxHydHz;                          col[6].i = i;   col[6].j = j;   col[6].k = k+1;
252:           MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
253:         }
254:       }
255:     }
256:   }
257:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
258:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
259:   if (user->bcType == NEUMANN) {
260:     MatNullSpace nullspace;

262:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
263:     MatSetNullSpace(jac,nullspace);
264:     MatNullSpaceDestroy(&nullspace);
265:   }
266:   return(0);
267: }