Actual source code: ex34.c

petsc-master 2018-01-17
Report Typos and Errors
  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;
 38:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,d,dof;
 39:   PetscScalar    Hx,Hy,Hz;
 40:   PetscScalar    ****array;
 41:   Vec            x,b,r;
 42:   Mat            J;

 44:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 45:   dof  = 1;
 46:   PetscOptionsGetInt(NULL,NULL,"-da_dof",&dof,NULL);
 47:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 48:   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);
 49:   DMSetFromOptions(da);
 50:   DMSetUp(da);
 51:   DMDASetInterpolationType(da, DMDA_Q0);

 53:   KSPSetDM(ksp,da);

 55:   KSPSetComputeRHS(ksp,ComputeRHS,NULL);
 56:   KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
 57:   KSPSetFromOptions(ksp);
 58:   KSPSolve(ksp,NULL,NULL);
 59:   KSPGetSolution(ksp,&x);
 60:   KSPGetRhs(ksp,&b);
 61:   KSPGetOperators(ksp,NULL,&J);
 62:   VecDuplicate(b,&r);

 64:   MatMult(J,x,r);
 65:   VecAXPY(r,-1.0,b);
 66:   VecNorm(r,NORM_2,&norm);
 67:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

 69:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
 70:   Hx   = 1.0 / (PetscReal)(mx);
 71:   Hy   = 1.0 / (PetscReal)(my);
 72:   Hz   = 1.0 / (PetscReal)(mz);
 73:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
 74:   DMDAVecGetArrayDOF(da, x, &array);

 76:   for (k=zs; k<zs+zm; k++) {
 77:     for (j=ys; j<ys+ym; j++) {
 78:       for (i=xs; i<xs+xm; i++) {
 79:         for (d=0; d<dof; d++) {
 80:           array[k][j][i][d] -=
 81:             PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
 82:             PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
 83:             PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
 84:         }
 85:       }
 86:     }
 87:   }
 88:   DMDAVecRestoreArrayDOF(da, x, &array);
 89:   VecAssemblyBegin(x);
 90:   VecAssemblyEnd(x);

 92:   VecNorm(x,NORM_INFINITY,&norm);
 93:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)norm);
 94:   VecNorm(x,NORM_1,&norm);
 95:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
 96:   VecNorm(x,NORM_2,&norm);
 97:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));

 99:   VecDestroy(&r);
100:   KSPDestroy(&ksp);
101:   DMDestroy(&da);
102:   PetscFinalize();
103:   return ierr;
104: }

106: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
107: {
109:   PetscInt       d,dof,i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
110:   PetscScalar    Hx,Hy,Hz;
111:   PetscScalar    ****array;
112:   DM             da;
113:   MatNullSpace   nullspace;

116:   KSPGetDM(ksp,&da);
117:   DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,&dof,0,0,0,0,0);
118:   Hx   = 1.0 / (PetscReal)(mx);
119:   Hy   = 1.0 / (PetscReal)(my);
120:   Hz   = 1.0 / (PetscReal)(mz);
121:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
122:   DMDAVecGetArrayDOF(da, b, &array);
123:   for (k=zs; k<zs+zm; k++) {
124:     for (j=ys; j<ys+ym; j++) {
125:       for (i=xs; i<xs+xm; i++) {
126:         for (d=0; d<dof; d++) {
127:           array[k][j][i][d] = 12 * PETSC_PI * PETSC_PI
128:                            * PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
129:                            * PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
130:                            * PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
131:                            * Hx * Hy * Hz;
132:         }
133:       }
134:     }
135:   }
136:   DMDAVecRestoreArrayDOF(da, b, &array);
137:   VecAssemblyBegin(b);
138:   VecAssemblyEnd(b);

140:   /* force right hand side to be consistent for singular matrix */
141:   /* note this is really a hack, normally the model would provide you with a consistent right handside */

143:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
144:   MatNullSpaceRemove(nullspace,b);
145:   MatNullSpaceDestroy(&nullspace);
146:   return(0);
147: }


150: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac, void *ctx)
151: {
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;

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,JJ2;

249:     MatComputeExplicitOperator(jac,&JJ);
250:     MatConvert(JJ,MATAIJ,MAT_INITIAL_MATRIX,&JJ2);
251:     MatChop(JJ2,1.e-8);
252:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_MATLAB);
253:     MatView(JJ2,PETSC_VIEWER_STDOUT_WORLD);
254:     MatDestroy(&JJ2);
255:     MatDestroy(&JJ);
256:   }
257:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
258:   MatSetNullSpace(J,nullspace);
259:   MatNullSpaceDestroy(&nullspace);
260:   return(0);
261: }