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