Actual source code: ex50.c

  1: /*   DMMG/KSP solving a system of linear equations.
  2:      Poisson equation in 2D: 

  4:      div(grad p) = f,  0 < x,y < 1
  5:      with
  6:        forcing function f = -cos(m*pi*x)*cos(n*pi*y),
  7:        Neuman boundary conditions
  8:         dp/dx = 0 for x = 0, x = 1.
  9:         dp/dy = 0 for y = 0, y = 1.

 11:      Contributed by Michael Boghosian <boghmic@iit.edu>, 2008,
 12:          based on petsc/src/ksp/ksp/examples/tutorials/ex29.c and ex32.c

 14:      Example of Usage: 
 15:           ./ex50 -mglevels 3 -ksp_monitor -M 3 -N 3 -ksp_view -da_view_draw -draw_pause -1 
 16:           ./ex50 -M 100 -N 100 -mglevels 1 -mg_levels_0_pc_factor_levels <ilu_levels> -ksp_monitor -cmp_solu
 17:           ./ex50 -M 100 -N 100 -mglevels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor -cmp_solu
 18:           mpiexec -n 4 ./ex50 -M 3 -N 3 -ksp_monitor -ksp_view -mglevels 10 -log_summary
 19: */

 21: static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";

 23: #include <petscdmda.h>
 24: #include <petscksp.h>
 25: #include <petscpcmg.h>
 26: #include <petscdmmg.h>
 27: #include <petscsys.h>
 28: #include <petscvec.h>


 35: typedef enum {DIRICHLET, NEUMANN} BCType;

 37: typedef struct {
 38:   PetscScalar  uu, tt;
 39:   BCType       bcType;
 40: } UserContext;

 44: int main(int argc,char **argv)
 45: {
 46:   DMMG           *dmmg;
 47:   DM             da;
 48:   UserContext    user;
 49:   PetscInt       l, bc, mglevels, M, N, stages[3];
 50:   PetscReal      norm;
 52:   PetscMPIInt    rank,nproc;
 53:   PetscBool      flg;
 54: 
 55:   PetscInitialize(&argc,&argv,(char *)0,help);
 56:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 57:   MPI_Comm_rank(PETSC_COMM_WORLD,&nproc);
 58:   PetscLogStageRegister("DMMG Setup",&stages[0]);
 59:   PetscLogStageRegister("DMMG Solve",&stages[1]);
 60: 
 61:   PetscLogStagePush(stages[0]);   /* Start DMMG Setup */
 62: 
 63:   /* SET VARIABLES: */
 64:   mglevels = 1;
 65:   PetscOptionsGetInt(PETSC_NULL,"-mglevels",&mglevels,PETSC_NULL);
 66:   M = 11;        /* number of grid points in x dir. on coarse grid */
 67:   N = 11;        /* number of grid points in y dir. on coarse grid */
 68:   PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
 69:   PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
 70: 
 71:   DMMGCreate(PETSC_COMM_WORLD,mglevels,PETSC_NULL,&dmmg);
 72:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);
 73:   DMMGSetDM(dmmg,(DM)da);
 74:   DMDestroy(&da);
 75: 
 76:   /* Set user contex */
 77:   user.uu = 1.0;
 78:   user.tt = 1.0;
 79:   bc   = (PetscInt)NEUMANN; // Use Neumann Boundary Conditions
 80:   user.bcType = (BCType)bc;
 81:   for (l = 0; l < DMMGGetLevels(dmmg); l++) {
 82:     DMMGSetUser(dmmg,l,&user);
 83:   }
 84: 
 85:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
 86:   if (user.bcType == NEUMANN){
 87:      DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
 88:   }
 89:   PetscLogStagePop();   /* Finish DMMG Setup */

 91:   /* DMMG SOLVE: */
 92:   PetscLogStagePush(stages[1]);   /* Start DMMG Solve */
 93:   DMMGSolve(dmmg);
 94:   PetscLogStagePop();   /* Finish DMMG Solve */

 96:   /* Compare solution with the true solution p */
 97:   PetscOptionsHasName(PETSC_NULL, "-cmp_solu", &flg);
 98:   if (flg){
 99:     Vec   x,p;
100:     if (mglevels != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"mglevels must equls 1 for comparison");
101:     if (nproc > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"num of proc must equls 1 for comparison");
102:     x = DMMGGetx(dmmg);
103:     VecDuplicate(x,&p);
104:     ComputeTrueSolution(dmmg,p);

106:     VecAXPY(p,-1.0,x);  /* p <- (p-x) */
107:     VecNorm(p, NORM_2, &norm);
108:     if (!rank){printf("| solu_compt - solu_true | = %g\n",norm);}
109:     VecDestroy(&p);
110:   }
111: 
112:   DMMGDestroy(dmmg);
113:   PetscFinalize();
114:   return 0;
115: }

117: // COMPUTE RHS:--------------------------------------------------------------
120: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
121: {
122:   DM             da = dmmg->dm;
123:   UserContext    *user = (UserContext *) dmmg->user;
125:   PetscInt       i, j, M, N, xm ,ym ,xs, ys;
126:   PetscScalar    Hx, Hy, pi, uu, tt;
127:   PetscScalar    **array;

130:   DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
131:   uu = user->uu; tt = user->tt;
132:   pi = 4*atan(1.0);
133:   Hx   = 1.0/(PetscReal)(M);
134:   Hy   = 1.0/(PetscReal)(N);

136:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);  // Fine grid
137:   //printf(" M N: %d %d; xm ym: %d %d; xs ys: %d %d\n",M,N,xm,ym,xs,ys);
138:   DMDAVecGetArray(da, b, &array);
139:   for (j=ys; j<ys+ym; j++){
140:     for(i=xs; i<xs+xm; i++){
141:       array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*cos(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
142:     }
143:   }
144:   DMDAVecRestoreArray(da, b, &array);
145:   VecAssemblyBegin(b);
146:   VecAssemblyEnd(b);

148:   /* force right hand side to be consistent for singular matrix */
149:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
150:   if (user->bcType == NEUMANN) {
151:     MatNullSpace nullspace;
152:     KSPGetNullSpace(dmmg->ksp,&nullspace);
153:     MatNullSpaceRemove(nullspace,b,PETSC_NULL);
154:   }
155:   return(0);
156: }

158: // COMPUTE JACOBIAN:--------------------------------------------------------------
161: PetscErrorCode ComputeJacobian(DMMG dmmg, Mat J, Mat jac)
162: {
163:   DM             da =  dmmg->dm;
164:   UserContext    *user = (UserContext *) dmmg->user;
166:   PetscInt       i, j, M, N, xm, ym, xs, ys, num, numi, numj;
167:   PetscScalar    v[5], Hx, Hy, HydHx, HxdHy;
168:   MatStencil     row, col[5];

171:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
172:   Hx    = 1.0 / (PetscReal)(M);
173:   Hy    = 1.0 / (PetscReal)(N);
174:   HxdHy = Hx/Hy;
175:   HydHx = Hy/Hx;
176:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
177:   for (j=ys; j<ys+ym; j++){
178:     for(i=xs; i<xs+xm; i++){
179:       row.i = i; row.j = j;
180: 
181:       if (i==0 || j==0 || i==M-1 || j==N-1) {
182:         if (user->bcType == DIRICHLET){
183:           SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
184:         } else if (user->bcType == NEUMANN){
185:           num=0; numi=0; numj=0;
186:           if (j!=0) {
187:             v[num] = -HxdHy;              col[num].i = i;   col[num].j = j-1;
188:             num++; numj++;
189:           }
190:           if (i!=0) {
191:             v[num] = -HydHx;              col[num].i = i-1; col[num].j = j;
192:             num++; numi++;
193:           }
194:           if (i!=M-1) {
195:             v[num] = -HydHx;              col[num].i = i+1; col[num].j = j;
196:             num++; numi++;
197:           }
198:           if (j!=N-1) {
199:             v[num] = -HxdHy;              col[num].i = i;   col[num].j = j+1;
200:             num++; numj++;
201:           }
202:           v[num] = ( (PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx ); col[num].i = i;   col[num].j = j;
203:           num++;
204:           MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
205:         }
206:       } else {
207:         v[0] = -HxdHy;              col[0].i = i;   col[0].j = j-1;
208:         v[1] = -HydHx;              col[1].i = i-1; col[1].j = j;
209:         v[2] = 2.0*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
210:         v[3] = -HydHx;              col[3].i = i+1; col[3].j = j;
211:         v[4] = -HxdHy;              col[4].i = i;   col[4].j = j+1;
212:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
213:       }
214:     }
215:   }
216:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
217:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
218:   return(0);
219: }

221: // COMPUTE TrueSolution:--------------------------------------------------------------
224: PetscErrorCode ComputeTrueSolution(DMMG *dmmg, Vec b)
225: {
226:   DM             da = (*dmmg)->dm;
227:   UserContext    *user = (UserContext *) (*dmmg)->user;
229:   PetscInt       i, j, M, N, xm ,ym ,xs, ys;
230:   PetscScalar    Hx, Hy, pi, uu, tt, cc;
231:   PetscScalar    *array;

234:   DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);  /* level_0 ! */
235:   //printf("ComputeTrueSolution - M N: %d %d;\n",M,N);

237:   uu = user->uu; tt = user->tt;
238:   pi = 4*atan(1.0);
239:   cc = -1.0/( (uu*pi)*(uu*pi) + (tt*pi)*(tt*pi) );
240:   Hx   = 1.0/(PetscReal)(M);
241:   Hy   = 1.0/(PetscReal)(N);

243:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
244:   VecGetArray(b, &array);
245:   for (j=ys; j<ys+ym; j++){
246:     for(i=xs; i<xs+xm; i++){
247:       array[j*xm+i] = cos(uu*pi*((PetscReal)i+0.5)*Hx)*cos(tt*pi*((PetscReal)j+0.5)*Hy)*cc;
248:     }
249:   }
250:   VecRestoreArray(b, &array);
251:   return(0);
252: }

254: // VECVIEW_VTK:--------------------------------------------------------------
257: PetscErrorCode VecView_VTK(Vec x, const char filename[], const char bcName[])
258: {
259:   MPI_Comm           comm;
260:   DM                 da;
261:   Vec                coords;
262:   PetscViewer        viewer;
263:   PetscScalar        *array, *values;
264:   PetscInt           nn, NN, maxn, M, N;
265:   PetscInt           i, p, dof = 1;
266:   MPI_Status         status;
267:   PetscMPIInt        rank, size, tag;
268:   PetscErrorCode     ierr;
269:   PetscBool          flg;

272:   PetscObjectGetComm((PetscObject) x, &comm);
273:   PetscViewerASCIIOpen(comm, filename, &viewer);

275:   VecGetSize(x, &NN);
276:   VecGetLocalSize(x, &nn);
277:   PetscObjectQuery((PetscObject) x, "DM", (PetscObject *) &da);
278:   if (!da) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
279:   PetscTypeCompare((PetscObject)da,DMDA,&flg);
280:   if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");

282:   DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,&dof,0,0,0,0,0);

284:   PetscViewerASCIIPrintf(viewer, "# vtk DataFile Version 2.0\n");
285:   PetscViewerASCIIPrintf(viewer, "Inhomogeneous Poisson Equation with %s boundary conditions\n", bcName);
286:   PetscViewerASCIIPrintf(viewer, "ASCII\n");
287:   // get coordinates of nodes
288:   DMDAGetCoordinates(da, &coords);
289:   if (!coords) {
290:     DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
291:     DMDAGetCoordinates(da, &coords);
292:   }
293:   PetscViewerASCIIPrintf(viewer, "DATASET RECTILINEAR_GRID\n");
294:   PetscViewerASCIIPrintf(viewer, "DIMENSIONS %d %d %d\n", M, N, 1);
295:   VecGetArray(coords, &array);
296:   PetscViewerASCIIPrintf(viewer, "X_COORDINATES %d double\n", M);
297:   for(i = 0; i < M; i++) {
298:     PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2]));
299:   }
300:   PetscViewerASCIIPrintf(viewer, "\n");
301:   PetscViewerASCIIPrintf(viewer, "Y_COORDINATES %d double\n", N);
302:   for(i = 0; i < N; i++) {
303:     PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*M*2+1]));
304:   }
305:   PetscViewerASCIIPrintf(viewer, "\n");
306:   PetscViewerASCIIPrintf(viewer, "Z_COORDINATES %d double\n", 1);
307:   PetscViewerASCIIPrintf(viewer, "%G\n", 0.0);
308:   VecRestoreArray(coords, &array);
309:   PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", N);
310:   PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", 1);
311:   PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
312:   VecGetArray(x, &array);

314:   // Determine maximum message to arrive:
315:   MPI_Comm_rank(comm, &rank);
316:   MPI_Comm_size(comm, &size) ;
317:   MPI_Reduce(&nn, &maxn, 1, MPIU_INT, MPI_MAX, 0, comm);
318:   tag  = ((PetscObject) viewer)->tag;
319:   if (!rank) {
320:     PetscMalloc((maxn+1) * sizeof(PetscScalar), &values);
321:     for(i = 0; i < nn; i++) {
322:       PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
323:     }
324:     for(p = 1; p < size; p++) {
325:       MPI_Recv(values, (PetscMPIInt) nn, MPIU_SCALAR, p, tag, comm, &status);
326:       MPI_Get_count(&status, MPIU_SCALAR, &nn);
327:       for(i = 0; i < nn; i++) {
328:         PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
329:       }
330:     }
331:     PetscFree(values);
332:   } else {
333:     MPI_Send(array, nn, MPIU_SCALAR, 0, tag, comm);
334:   }
335:   VecRestoreArray(x, &array);
336:   PetscViewerFlush(viewer);
337:   PetscViewerDestroy(&viewer);
338: 
339:   return(0);
340:   }