Actual source code: ex47.c

  2: /*
  3: Laplacian in 3D. Modeled by the partial differential equation

  5:    - Laplacian u = 1,0 < x,y,z < 1,

  7: with boundary conditions

  9:    u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.

 11:    This uses multigrid to solve the linear system

 13: */

 15: static char help[] = "Solves 3D Laplacian using multigrid.\n\n";

 17: #include <petscdmda.h>
 18: #include <petscksp.h>
 19: #include <petscdmmg.h>


 27: PetscReal L[3] = {1.0, 1.0, 1.0};

 31: int main(int argc,char **argv)
 32: {
 34:   DMMG           *dmmg;
 35:   PetscReal      norm, normTotal;
 36:   DM             da;
 37:   Vec            phi, phiRhs;

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

 41:   DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
 42:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-3,-3,-3,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
 43:   DMMGSetDM(dmmg,(DM)da);
 44:   DMDestroy(&da);

 46:   DMMGSetKSP(dmmg,(PetscErrorCode (*)(DMMG, Vec)) ComputeRHS,ComputeMatrix);
 47:   //DMMGSetNullSpace(dmmg, PETSC_TRUE, 0, PETSC_NULL);

 49:   DMMGSetUp(dmmg);
 50:   DMMGSolve(dmmg);

 52:   VecDuplicate(DMMGGetx(dmmg), &phi);
 53:   VecDuplicate(DMMGGetx(dmmg), &phiRhs);
 54:   ComputeRHS(dmmg[0], phiRhs, PETSC_TRUE);
 55:   Solve_FFT(DMMGGetDM(dmmg), phiRhs, phi);

 57:   Vec       stddev;
 58:   PetscReal s;
 59:   PetscInt  p;

 61:   CalculateXYStdDev(da, DMMGGetRHS(dmmg), &stddev);
 62:   VecMax(stddev, &p, &s);
 63:   if (s > 1.0e-10) {
 64:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "RHS Homogeneity violation, std deviation %g z %d", s, p);
 65:   } else {
 66:     PetscPrintf(PETSC_COMM_WORLD, "FD RHS Homogeneity verification\n");
 67:     PetscPrintf(PETSC_COMM_WORLD, "    std deviation   %g\n", s);
 68:   }
 69:   VecDestroy(&stddev);
 70:   CalculateXYStdDev(da, DMMGGetx(dmmg), &stddev);
 71:   VecMax(stddev, &p, &s);
 72:   if (s > 1.0e-5) {
 73:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "FD Homogeneity violation, std deviation %g z %d", s, p);
 74:   } else {
 75:     PetscPrintf(PETSC_COMM_WORLD, "FD Solution Homogeneity verification\n");
 76:     PetscPrintf(PETSC_COMM_WORLD, "    std deviation   %g\n", s);
 77:   }
 78:   VecDestroy(&stddev);
 79:   CalculateXYStdDev(da, phi, &stddev);
 80:   VecMax(stddev, &p, &s);
 81:   if (s > 1.0e-10) {
 82:     SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB, "FFT Homogeneity violation, std deviation %g z %d", s, p);
 83:   } else {
 84:     PetscPrintf(PETSC_COMM_WORLD, "FFT Solution Homogeneity verification\n");
 85:     PetscPrintf(PETSC_COMM_WORLD, "    std deviation   %g\n", s);
 86:   }
 87:   VecDestroy(&stddev);

 89:   PetscInt N;
 90:   Vec      tmp;
 91:   VecGetSize(phi, &N);
 92:   VecCreate(PETSC_COMM_WORLD, &tmp);
 93:   VecSetSizes(tmp, PETSC_DECIDE, N);
 94:   VecSetFromOptions(tmp);

 96:   VecCopy(DMMGGetx(dmmg), tmp);
 97:   VecView(tmp, PETSC_VIEWER_DRAW_WORLD);
 98:   VecViewCenterSingle(da, DMMGGetx(dmmg), PETSC_VIEWER_DRAW_WORLD, "FD Solution", -1, -1);
 99:   VecViewCenterSingle(da, DMMGGetx(dmmg), PETSC_VIEWER_STDOUT_WORLD, "FD Solution", -1, -1);
100:   MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
101:   VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
102:   VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
103:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);

105:   VecCopy(phi, tmp);
106:   VecView(tmp, PETSC_VIEWER_DRAW_WORLD);
107:   VecViewCenterSingle(da, phi, PETSC_VIEWER_DRAW_WORLD, "FFT Solution", -1, -1);
108:   VecViewCenterSingle(da, phi, PETSC_VIEWER_STDOUT_WORLD, "FFT Solution", -1, -1);
109:   VecAXPY(phi,-1.0,DMMGGetx(dmmg));
110:   VecNorm(phi,NORM_2,&norm);
111:   VecNorm(DMMGGetx(dmmg),NORM_2,&normTotal);
112:   PetscPrintf(PETSC_COMM_WORLD,"Error norm (FFT vs. FD) %G %G\n",norm, norm/normTotal);

114:   VecCopy(phi, tmp);
115:   VecView(tmp, PETSC_VIEWER_DRAW_WORLD);
116:   VecViewCenterSingle(da, phi, PETSC_VIEWER_DRAW_WORLD, "Error", -1, -1);
117:   VecViewCenterSingle(da, phi, PETSC_VIEWER_STDOUT_WORLD, "Error", -1, -1);
118:   VecDestroy(&tmp);

120:   DMMGDestroy(dmmg);
121:   PetscFinalize();
122:   return 0;
123: }

127: PetscErrorCode ComputeRHS(DMMG dmmg,Vec b, PetscBool  withBC= PETSC_TRUE)
128: {
129:   DM             da =  dmmg->dm;
130:   PetscInt       bathIndex;
131:   PetscScalar ***a;
132:   PetscScalar    sc;
133:   PetscInt       mx, my, mz, xm, ym, zm, xs, ys, zs, wallPos, i, j, k;

137:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
138:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
139:   sc   = 1.0/((mx-1)*(my-1)*(mz-1));
140:   //wallPos = (mz-1)/20;
141:   wallPos = -1;
142:   bathIndex = mz/2;
143:   DMDAVecGetArray(da, b, &a);
144:   for(k = zs; k < zs+zm; ++k) {
145:     for(j = ys; j < ys+ym; ++j) {
146:       for(i = xs; i < xs+xm; ++i) {
147:         if (k == bathIndex && withBC) {
148:           a[k][j][i] = 0.0;
149:         } else {
150:           if (k > wallPos) {
151:             a[k][j][i] = sc;
152:           } else {
153:             a[k][j][i] = 0.0;
154:           }
155:         }
156:       }
157:     }
158:   }
159:   DMDAVecRestoreArray(da, b, &a);
160:   return(0);
161: }
162: 
165: PetscErrorCode ComputeMatrix(DMMG dmmg,Mat jac,Mat B)
166: {
167:   DM             da = dmmg->dm;
168:   PetscInt       bathIndex;
170:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
171:   PetscScalar    v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
172:   MatStencil     row,col[7];

174:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
175:   Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
176:   //Hx = L[0] / (PetscReal)(mx-1); Hy = L[1] / (PetscReal)(my-1); Hz = L[2] / (PetscReal)(mz-1);
177:   //HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
178:   HxHydHz = 1.0/PetscSqr(Hz); HxHzdHy = 1.0/PetscSqr(Hy); HyHzdHx = 1.0/PetscSqr(Hx);
179:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
180:   bathIndex = mz/2;
182:   for (k=zs; k<zs+zm; k++){
183:     for (j=ys; j<ys+ym; j++){
184:       for(i=xs; i<xs+xm; i++){
185:         row.i = i; row.j = j; row.k = k;
186:         if (k == bathIndex) {
187:           v[0] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
188:           MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
189:         } else {
190:           v[0] = -HxHydHz;col[0].i = i; col[0].j = j; col[0].k = k-1;
191:           v[1] = -HxHzdHy;col[1].i = i; col[1].j = j-1; col[1].k = k;
192:           v[2] = -HyHzdHx;col[2].i = i-1; col[2].j = j; col[2].k = k;
193:           v[3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);col[3].i = row.i; col[3].j = row.j; col[3].k = row.k;
194:           v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
195:           v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
196:           v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
197:           MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);
198:         }
199:       }
200:     }
201:   }
202:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
203:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
204:   return(0);
205: }


210: PetscErrorCode Solve_FFT(DM da, Vec rhs, Vec phi)
211: {
212:   PetscReal      h[3];
213:   PetscInt       dim[3];
214:   PetscReal      scale, sc;
215:   PetscScalar ***rhsHatArray;
216:   PetscScalar ***phiHatArray;
217:   Mat            F;
218:   Vec            rhsHat, phiHat;
219:   PetscInt       M, N, P, xm, ym, zm, xs, ys, zs;
221: 
223:   DMDAGetInfo(da, 0, &M, &N, &P, 0, 0, 0, 0,0,0,0,0,0);
224:   DMDAGetInfo(da, 0, &dim[2], &dim[1], &dim[0], 0, 0, 0, 0,0,0,0,0,0);
225:   MatCreateSeqFFTW(PETSC_COMM_WORLD, 3, dim, &F);
226:   DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm);
227:   h[0] = L[0]/(M - 1);
228:   h[1] = L[1]/(N - 1);
229:   h[2] = L[2]/(P - 1);
230:   scale = 1.0/((PetscReal) M*N*P);
231:   sc    = (M-1)*(N-1)*(P-1);
232:   DMGetGlobalVector(da, &rhsHat);
233:   MatMult(F, rhs, rhsHat);
234:   DMGetGlobalVector(da, &phiHat);
235:   DMDAVecGetArray(da, rhsHat, &rhsHatArray);
236:   DMDAVecGetArray(da, phiHat, &phiHatArray);
237:   for(PetscInt k = zs; k < zs+zm; ++k) {
238:     PetscReal kz = k <= P/2 ? 2.0*M_PI*k/(P) : -2.0*M_PI*(P-k)/(P);

240:     for(PetscInt j = ys; j < ys+ym; ++j) {
241:       PetscReal ky = j <= N/2 ? 2.0*M_PI*j/(N) : -2.0*M_PI*(N-j)/(N);

243:       for(PetscInt i = xs; i < xs+xm; ++i) {
244:         PetscReal kx = i <= M/2 ? 2.0*M_PI*i/(M) : -2.0*M_PI*(M-i)/(M);
245:         PetscScalar charge = 0.0;

247: #if 0
248:         for(PetscInt sp = 0; sp < geomOptions->numSpecies; ++sp) {
249:           charge += (e*e/epsilon)*z[sp]*rhoHatArray[k][j][i].v[sp];
250:         }
251: #else
252:         charge = rhsHatArray[k][j][i];
253: #endif
254:         PetscReal denom = 2.0*((1.0-cos(kx))/PetscSqr(h[0]) + (1.0-cos(ky))/PetscSqr(h[1]) + (1.0-cos(kz))/PetscSqr(h[2]));
255:         // Omit the zeroth moment
256:         if (kx == 0 && ky == 0 && kz == 0) {
257:           phiHatArray[k][j][i] = 0.0;
258:         } else {
259:           phiHatArray[k][j][i] = charge/denom; // Changed units by scaling by e
260:         }
261:         if (PetscIsInfOrNanScalar(phiHatArray[k][j][i])) {
262:           SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FP, "Nan or inf at phiHat[%d][%d][%d]: %g ", k, j, i, PetscRealPart(phiHatArray[k][j][i]));
263:         }
264:       }
265:     }
266:   }
267:   DMDAVecRestoreArray(da, phiHat, &phiHatArray);
268:   MatMultTranspose(F, phiHat, phi);
269:   VecScale(phi, scale);
270:   DMRestoreGlobalVector(da, &phiHat);
271:   DMRestoreGlobalVector(da, &rhsHat);
272:   MatDestroy(&F);

274:   // Force potential in the bath to be 0
275:   PetscInt       bathIndex[3] = {0, 0, P/2};
276:   PetscScalar ***phiArray;
277:   PetscScalar    bathPotential;

279:   DMDAVecGetArray(da, phi, &phiArray);
280:   bathPotential = phiArray[bathIndex[2]][bathIndex[1]][bathIndex[0]];
281:   DMDAVecRestoreArray(da, phi, &phiArray);
282:   VecShift(phi, -bathPotential);
283:   return(0);
284: }

288: PetscErrorCode CalculateXYStdDev(DM da, Vec v, Vec *std) {
289:   DMDALocalInfo    info;
290:   MPI_Comm       comm;
291:   PetscScalar ***a;
292:   PetscScalar   *r;
294: 
296:   PetscObjectGetComm((PetscObject) da, &comm);
297:   DMDAGetLocalInfo(da, &info);
298:   DMDAVecGetArray(da, v, &a);
299:     VecCreate(comm, std);
300:     VecSetSizes(*std, info.zm - info.zs, info.mz);
301:     VecSetFromOptions(*std);
302:     VecGetArray(*std, &r);
303:     for(PetscInt k = 0; k < info.mz; ++k) {
304:       PetscScalar avg = 0.0, var = 0.0;

306:       for(PetscInt j = 0; j < info.my; ++j) {
307:         for(PetscInt i = 0; i < info.mx; ++i) {
308:           avg += a[k][j][i];
309:         }
310:       }
311:       avg /= (info.mx*info.my);
312:       for(PetscInt j = 0; j < info.my; ++j) {
313:         for(PetscInt i = 0; i < info.mx; ++i) {
314:           var += PetscSqr(a[k][j][i] - avg);
315:         }
316:       }
317:       r[k] = PetscSqrtScalar(var);
318:         }
319:     VecRestoreArray(*std, &r);
320:   DMDAVecRestoreArray(da, v, &a);
321:   return(0);
322: }

326: PetscErrorCode VecViewCenterSingle(DM da, Vec v, PetscViewer viewer, const char name[], PetscInt i, PetscInt j)
327: {
328:   DMDALocalInfo    info;
329:   MPI_Comm       comm;
330:   Vec            c;
331:   PetscScalar ***a;
332:   PetscScalar   *b;

336:   PetscObjectGetComm((PetscObject) da, &comm);
337:   if (i < 0) {
338:     PetscPrintf(comm, "Viewing %s\n", name);
339:   } else {
340:     if (j < 0) {
341:       PetscPrintf(comm, "Viewing %s[%d]\n", name, i);
342:     } else {
343:       PetscPrintf(comm, "Viewing %s[%d,%d]\n", name, i, j);
344:     }
345:   }
346:   DMDAGetLocalInfo(da, &info);
347:   VecCreate(comm, &c);
348:   VecSetSizes(c, info.zm - info.zs, info.mz);
349:   VecSetFromOptions(c);
350:   DMDAVecGetArray(da, v, &a);
351:   VecGetArray(c, &b);
352:   for(PetscInt k = 0, i = info.mx/2, j = info.my/2; k < info.mz; ++k) {
353:     b[k] = a[k][j][i];
354:   }
355:   DMDAVecRestoreArray(da, v, &a);
356:   VecRestoreArray(c, &b);
357:   VecView(c, viewer);
358:   VecDestroy(&c);
359:   return(0);
360: }