Actual source code: ex34.c

petsc-master 2016-06-26
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*);

 34: int main(int argc,char **argv)
 35: {
 36:   KSP            ksp;
 37:   DM             da;
 38:   PetscReal      norm;
 40:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
 41:   PetscScalar    Hx,Hy,Hz;
 42:   PetscScalar    ***array;
 43:   Vec            x,b,r;
 44:   Mat            J;

 46:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 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,1,1,0,0,0,&da);
 49:   DMDASetInterpolationType(da, DMDA_Q0);

 51:   KSPSetDM(ksp,da);

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

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

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

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

 88:   VecNorm(x,NORM_INFINITY,&norm);
 89:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)norm);
 90:   VecNorm(x,NORM_1,&norm);
 91:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
 92:   VecNorm(x,NORM_2,&norm);
 93:   PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));

 95:   VecDestroy(&r);
 96:   KSPDestroy(&ksp);
 97:   DMDestroy(&da);
 98:   PetscFinalize();
 99:   return ierr;
100: }

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

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

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

139:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
140:   MatNullSpaceRemove(nullspace,b);
141:   MatNullSpaceDestroy(&nullspace);
142:   return(0);
143: }


148: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac, void *ctx)
149: {
151:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
152:   PetscScalar    v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
153:   MatStencil     row, col[7];
154:   DM             da;
155:   MatNullSpace   nullspace;

158:   KSPGetDM(ksp,&da);
159:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
160:   Hx      = 1.0 / (PetscReal)(mx);
161:   Hy      = 1.0 / (PetscReal)(my);
162:   Hz      = 1.0 / (PetscReal)(mz);
163:   HyHzdHx = Hy*Hz/Hx;
164:   HxHzdHy = Hx*Hz/Hy;
165:   HxHydHz = Hx*Hy/Hz;
166:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
167:   for (k=zs; k<zs+zm; k++) {
168:     for (j=ys; j<ys+ym; j++) {
169:       for (i=xs; i<xs+xm; i++) {
170:         row.i = i; row.j = j; row.k = k;
171:         if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
172:           num = 0; numi=0; numj=0; numk=0;
173:           if (k!=0) {
174:             v[num]     = -HxHydHz;
175:             col[num].i = i;
176:             col[num].j = j;
177:             col[num].k = k-1;
178:             num++; numk++;
179:           }
180:           if (j!=0) {
181:             v[num]     = -HxHzdHy;
182:             col[num].i = i;
183:             col[num].j = j-1;
184:             col[num].k = k;
185:             num++; numj++;
186:             }
187:           if (i!=0) {
188:             v[num]     = -HyHzdHx;
189:             col[num].i = i-1;
190:             col[num].j = j;
191:             col[num].k = k;
192:             num++; numi++;
193:           }
194:           if (i!=mx-1) {
195:             v[num]     = -HyHzdHx;
196:             col[num].i = i+1;
197:             col[num].j = j;
198:             col[num].k = k;
199:             num++; numi++;
200:           }
201:           if (j!=my-1) {
202:             v[num]     = -HxHzdHy;
203:             col[num].i = i;
204:             col[num].j = j+1;
205:             col[num].k = k;
206:             num++; numj++;
207:           }
208:           if (k!=mz-1) {
209:             v[num]     = -HxHydHz;
210:             col[num].i = i;
211:             col[num].j = j;
212:             col[num].k = k+1;
213:             num++; numk++;
214:           }
215:           v[num]     = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
216:           col[num].i = i;   col[num].j = j;   col[num].k = k;
217:           num++;
218:           MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
219:         } else {
220:           v[0] = -HxHydHz;                          col[0].i = i;   col[0].j = j;   col[0].k = k-1;
221:           v[1] = -HxHzdHy;                          col[1].i = i;   col[1].j = j-1; col[1].k = k;
222:           v[2] = -HyHzdHx;                          col[2].i = i-1; col[2].j = j;   col[2].k = k;
223:           v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i;   col[3].j = j;   col[3].k = k;
224:           v[4] = -HyHzdHx;                          col[4].i = i+1; col[4].j = j;   col[4].k = k;
225:           v[5] = -HxHzdHy;                          col[5].i = i;   col[5].j = j+1; col[5].k = k;
226:           v[6] = -HxHydHz;                          col[6].i = i;   col[6].j = j;   col[6].k = k+1;
227:           MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
228:         }
229:       }
230:     }
231:   }
232:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
233:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
234:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
235:   MatSetNullSpace(J,nullspace);
236:   MatNullSpaceDestroy(&nullspace);
237:   return(0);
238: }