Actual source code: ex46.c

  1: static char help[] = "Surface processes in geophysics.\n\n";

  3: /*T
  4:    Concepts: SNES^parallel Surface process example
  5:    Concepts: DMDA^using distributed arrays;
  6:    Concepts: IS coloirng types;
  7:    Processors: n
  8: T*/

 10: /* ------------------------------------------------------------------------

 12:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 13:     the partial differential equation
 14:   
 15:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 16:   
 17:     with boundary conditions
 18:    
 19:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 20:   
 21:     A finite difference approximation with the usual 5-point stencil
 22:     is used to discretize the boundary value problem to obtain a nonlinear 
 23:     system of equations.

 25:     Program usage:  mpiexec -n <procs> ex5 [-help] [all PETSc options] 
 26:      e.g.,
 27:       ./ex5 -fd_jacobian -mat_fd_coloring_view_draw -draw_pause -1
 28:       mpiexec -n 2 ./ex5 -fd_jacobian_ghosted -log_summary

 30:   ------------------------------------------------------------------------- */

 32: /* 
 33:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 34:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 35:    file automatically includes:
 36:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 37:      petscmat.h - matrices
 38:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 39:      petscviewer.h - viewers               petscpc.h  - preconditioners
 40:      petscksp.h   - linear solvers
 41: */
 42: #include <petscdmmg.h>
 43: #include <petscsnes.h>

 45: /* 
 46:    User-defined application context - contains data needed by the 
 47:    application-provided call-back routines, FormJacobianLocal() and
 48:    FormFunctionLocal().
 49: */
 50: typedef struct {
 51:   DM          da; /* distributed array data structure */
 52:   PassiveReal D;  /* The diffusion coefficient */
 53:   PassiveReal K;  /* The advection coefficient */
 54:   PetscInt    m;  /* Exponent for A */
 55: } AppCtx;

 57: /* 
 58:    User-defined routines
 59: */

 66: int main(int argc,char **argv)
 67: {
 68:   DMMG                  *dmmg;
 69:   SNES                   snes;                 /* nonlinear solver */
 70:   AppCtx                 user;                 /* user-defined work context */
 71:   PetscInt               its;                  /* iterations for convergence */
 72:   PetscErrorCode         ierr;

 74:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75:      Initialize program
 76:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

 80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 81:      Initialize problem parameters
 82:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 83:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Surface Process Problem Options", "DMMG");
 84:     user.D = 1.0;
 85:     PetscOptionsReal("-D", "The diffusion coefficient D", __FILE__, user.D, &user.D, PETSC_NULL);
 86:     user.K = 1.0;
 87:     PetscOptionsReal("-K", "The advection coefficient K", __FILE__, user.K, &user.K, PETSC_NULL);
 88:     user.m = 1;
 89:     PetscOptionsInt("-m", "The exponent for A", __FILE__, user.m, &user.m, PETSC_NULL);
 90:   PetscOptionsEnd();

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:      Create distributed array (DMDA) to manage parallel grid and vectors
 94:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
 96:                     1,1,PETSC_NULL,PETSC_NULL,&user.da);
 97:   DMDASetUniformCoordinates(user.da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
 98:   DMMGCreate(PETSC_COMM_WORLD, 1, &user, &dmmg);
 99:   DMMGSetDM(dmmg, (DM) user.da);

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102:      Set local function evaluation routine
103:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104:   DMMGSetSNESLocal(dmmg, FormFunctionLocal, FormJacobianLocal, 0, 0);

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Customize solver; set runtime options
108:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109:   DMMGSetFromOptions(dmmg);

111:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:      Form initial guess
113:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114:   //FormInitialGuess(&user,DMMGGetx(dmmg));

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Solve nonlinear system
118:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119:   DMMGSolve(dmmg);
120:   snes = DMMGGetSNES(dmmg);
121:   SNESGetIterationNumber(snes,&its);

123:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);

127:   VecAssemblyBegin(DMMGGetx(dmmg));
128:   VecAssemblyEnd(DMMGGetx(dmmg));
129:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130:      Free work space.  All PETSc objects should be destroyed when they
131:      are no longer needed.
132:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

134:   DMMGDestroy(dmmg);
135:   DMDestroy(&user.da);

137:   PetscFinalize();
138:   return(0);
139: }
140: /* ------------------------------------------------------------------- */
143: /* 
144:    FormInitialGuess - Forms initial approximation.

146:    Input Parameters:
147:    user - user-defined application context
148:    X - vector

150:    Output Parameter:
151:    X - vector
152:  */
153: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
154: {
155:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
157:   PetscReal      D,temp1,temp,hx,hy;
158:   PetscScalar    **x;

161:   DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
162:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

164:   D      = user->D;
165:   hx     = 1.0/(PetscReal)(Mx-1);
166:   hy     = 1.0/(PetscReal)(My-1);
167:   temp1  = D/(D + 1.0);

169:   /*
170:      Get a pointer to vector data.
171:        - For default PETSc vectors, VecGetArray() returns a pointer to
172:          the data array.  Otherwise, the routine is implementation dependent.
173:        - You MUST call VecRestoreArray() when you no longer need access to
174:          the array.
175:   */
176:   DMDAVecGetArray(user->da,X,&x);

178:   /*
179:      Get local grid boundaries (for 2-dimensional DMDA):
180:        xs, ys   - starting grid indices (no ghost points)
181:        xm, ym   - widths of local grid (no ghost points)

183:   */
184:   DMDAGetCorners(user->da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

186:   /*
187:      Compute initial guess over the locally owned part of the grid
188:   */
189:   for (j=ys; j<ys+ym; j++) {
190:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
191:     for (i=xs; i<xs+xm; i++) {
192:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
193:         /* boundary conditions are all zero Dirichlet */
194:         x[j][i] = 0.0;
195:       } else {
196:         x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
197:       }
198:     }
199:   }

201:   /*
202:      Restore vector
203:   */
204:   DMDAVecRestoreArray(user->da,X,&x);

206:   return(0);
207: }

211: PetscScalar funcU(DMDACoor2d *coords)
212: {
213:   return coords->x + coords->y;
214: }

218: PetscScalar funcA(PetscScalar z, AppCtx *user)
219: {
220:   PetscScalar v = 1.0;
221:   PetscInt    i;

223:   for(i = 0; i < user->m; ++i) {
224:     v *= z;
225:   }
226:   return v;
227: }

231: PetscScalar funcADer(PetscScalar z, AppCtx *user)
232: {
233:   PetscScalar v = 1.0;
234:   PetscInt    i;

236:   for(i = 0; i < user->m-1; ++i) {
237:     v *= z;
238:   }
239:   return user->m*v;
240: }

244: /* 
245:    FormFunctionLocal - Evaluates nonlinear function, F(x).
246: */
247: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
248: {
249:   DM             coordDA;
250:   Vec            coordinates;
251:   DMDACoor2d     **coords;
252:   PetscScalar    u, ux, uy, uxx, uyy;
253:   PetscReal      D, K, hx, hy, hxdhy, hydhx;
254:   PetscInt       i,j;


259:   D      = user->D;
260:   K      = user->K;
261:   hx     = 1.0/(PetscReal)(info->mx-1);
262:   hy     = 1.0/(PetscReal)(info->my-1);
263:   hxdhy  = hx/hy;
264:   hydhx  = hy/hx;
265:   /*
266:      Compute function over the locally owned part of the grid
267:   */
268:   DMDAGetCoordinateDA(user->da, &coordDA);
269:   DMDAGetCoordinates(user->da, &coordinates);
270:   DMDAVecGetArray(coordDA, coordinates, &coords);
271:   for (j=info->ys; j<info->ys+info->ym; j++) {
272:     for (i=info->xs; i<info->xs+info->xm; i++) {
273:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
274:         f[j][i] = x[j][i];
275:       } else {
276:         u       = x[j][i];
277:         ux      = (x[j][i+1] - x[j][i])/hx;
278:         uy      = (x[j+1][i] - x[j][i])/hy;
279:         uxx     = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
280:         uyy     = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
281:         f[j][i] = D*(uxx + uyy) - (K*funcA(x[j][i], user)*sqrt(ux*ux + uy*uy) + funcU(&coords[j][i]))*hx*hy;
282:         if (PetscIsInfOrNanScalar(f[j][i])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", PetscRealPart(f[j][i]));
283:       }
284:     }
285:   }
286:   DMDAVecRestoreArray(coordDA, coordinates, &coords);
287:   PetscLogFlops(11*info->ym*info->xm);
288:   return(0);
289: }

293: /*
294:    FormJacobianLocal - Evaluates Jacobian matrix.
295: */
296: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat jac,AppCtx *user)
297: {
298:   MatStencil     col[5], row;
299:   PetscScalar    D, K, A, v[5], hx, hy, hxdhy, hydhx, ux, uy, normGradZ;
300:   PetscInt       i, j,k;

304:   D      = user->D;
305:   K      = user->K;
306:   hx     = 1.0/(PetscReal)(info->mx-1);
307:   hy     = 1.0/(PetscReal)(info->my-1);
308:   hxdhy  = hx/hy;
309:   hydhx  = hy/hx;

311:   /* 
312:      Compute entries for the locally owned part of the Jacobian.
313:       - Currently, all PETSc parallel matrix formats are partitioned by
314:         contiguous chunks of rows across the processors. 
315:       - Each processor needs to insert only elements that it owns
316:         locally (but any non-local elements will be sent to the
317:         appropriate processor during matrix assembly). 
318:       - Here, we set all entries for a particular row at once.
319:       - We can set matrix entries either using either
320:         MatSetValuesLocal() or MatSetValues(), as discussed above.
321:   */
322:   for (j=info->ys; j<info->ys+info->ym; j++) {
323:     for (i=info->xs; i<info->xs+info->xm; i++) {
324:       row.j = j; row.i = i;
325:       if (i == 0 || j == 0 || i == info->mx-1 || j == info->my-1) {
326:         /* boundary points */
327:         v[0] = 1.0;
328:         MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
329:       } else {
330:         /* interior grid points */
331:         ux        = (x[j][i+1] - x[j][i])/hx;
332:         uy        = (x[j+1][i] - x[j][i])/hy;
333:         normGradZ = sqrt(ux*ux + uy*uy);
334:         //PetscPrintf(PETSC_COMM_SELF, "i: %d j: %d normGradZ: %g\n", i, j, normGradZ);
335:         if (normGradZ < 1.0e-8) {
336:           normGradZ = 1.0e-8;
337:         }
338:         A         = funcA(x[j][i], user);

340:         v[0] = -D*hxdhy;                                                                          col[0].j = j - 1; col[0].i = i;
341:         v[1] = -D*hydhx;                                                                          col[1].j = j;     col[1].i = i-1;
342:         v[2] = D*2.0*(hydhx + hxdhy) + K*(funcADer(x[j][i], user)*normGradZ - A/normGradZ)*hx*hy; col[2].j = row.j; col[2].i = row.i;
343:         v[3] = -D*hydhx + K*A*hx*hy/(2.0*normGradZ);                                              col[3].j = j;     col[3].i = i+1;
344:         v[4] = -D*hxdhy + K*A*hx*hy/(2.0*normGradZ);                                              col[4].j = j + 1; col[4].i = i;
345:         for(k = 0; k < 5; ++k) {
346:           if (PetscIsInfOrNanScalar(v[k])) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FP, "Invalid residual: %g", PetscRealPart(v[k]));
347:         }
348:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
349:       }
350:     }
351:   }

353:   /* 
354:      Assemble matrix, using the 2-step process:
355:        MatAssemblyBegin(), MatAssemblyEnd().
356:   */
357:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
358:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
359:   /*
360:      Tell the matrix we will never add a new nonzero location to the
361:      matrix. If we do, it will generate an error.
362:   */
363:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
364:   return(0);
365: }