Actual source code: ex5.c

petsc-3.9.4 2018-09-11
Report Typos and Errors

  2: static char help[] = "Demonstrates Pattern Formation with Reaction-Diffusion Equations.\n";

  4: /*
  5:      Page 21, Pattern Formation with Reaction-Diffusion Equations

  7:         u_t = D1 (u_xx + u_yy)  - u*v^2 + gama(1 -u)
  8:         v_t = D2 (v_xx + v_yy)  + u*v^2 - (gamma + kappa)v

 10:     Unlike in the book this uses periodic boundary conditions instead of Neumann
 11:     (since they are easier for finite differences).
 12: */

 14: /*
 15:       Helpful runtime monitor options:
 16:            -ts_monitor_draw_solution
 17:            -draw_save -draw_save_movie

 19:       Helpful runtime linear solver options:
 20:            -pc_type mg -pc_mg_galerkin pmat -da_refine 1 -snes_monitor -ksp_monitor -ts_view  (note that these Jacobians are so well-conditioned multigrid may not be the best solver)

 22:       Point your browser to localhost:8080 to monitor the simulation
 23:            ./ex5  -ts_view_pre saws  -stack_view saws -draw_save -draw_save_single_file -x_virtual -ts_monitor_draw_solution -saws_root .

 25: */

 27: /*

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

 42: typedef struct {
 43:   PetscScalar u,v;
 44: } Field;

 46: typedef struct {
 47:   PetscReal D1,D2,gamma,kappa;
 48: } AppCtx;

 50: /*
 51:    User-defined routines
 52: */
 53: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*),InitialConditions(DM,Vec);
 54: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);

 56: int main(int argc,char **argv)
 57: {
 58:   TS             ts;                  /* ODE integrator */
 59:   Vec            x;                   /* solution */
 61:   DM             da;
 62:   AppCtx         appctx;

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:      Initialize program
 66:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 67:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 69:   appctx.D1    = 8.0e-5;
 70:   appctx.D2    = 4.0e-5;
 71:   appctx.gamma = .024;
 72:   appctx.kappa = .06;

 74:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 75:      Create distributed array (DMDA) to manage parallel grid and vectors
 76:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 77:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,65,65,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
 78:   DMSetFromOptions(da);
 79:   DMSetUp(da);
 80:   DMDASetFieldName(da,0,"u");
 81:   DMDASetFieldName(da,1,"v");

 83:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 84:      Extract global vectors from DMDA; then duplicate for remaining
 85:      vectors that are the same types
 86:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 87:   DMCreateGlobalVector(da,&x);

 89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 90:      Create timestepping solver context
 91:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 92:   TSCreate(PETSC_COMM_WORLD,&ts);
 93:   TSSetType(ts,TSARKIMEX);
 94:   TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);
 95:   TSSetDM(ts,da);
 96:   TSSetProblemType(ts,TS_NONLINEAR);
 97:   TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);
 98:   TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Set initial conditions
102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   InitialConditions(da,x);
104:   TSSetSolution(ts,x);

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Set solver options
108:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109:   TSSetMaxTime(ts,2000.0);
110:   TSSetTimeStep(ts,.0001);
111:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
112:   TSSetFromOptions(ts);

114:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
115:      Solve ODE system
116:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117:   TSSolve(ts,x);

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Free work space.  All PETSc objects should be destroyed when they
121:      are no longer needed.
122:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123:   VecDestroy(&x);
124:   TSDestroy(&ts);
125:   DMDestroy(&da);

127:   PetscFinalize();
128:   return ierr;
129: }
130: /* ------------------------------------------------------------------- */
131: /*
132:    RHSFunction - Evaluates nonlinear function, F(x).

134:    Input Parameters:
135: .  ts - the TS context
136: .  X - input vector
137: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

139:    Output Parameter:
140: .  F - function vector
141:  */
142: PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
143: {
144:   AppCtx         *appctx = (AppCtx*)ptr;
145:   DM             da;
147:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
148:   PetscReal      hx,hy,sx,sy;
149:   PetscScalar    uc,uxx,uyy,vc,vxx,vyy;
150:   Field          **u,**f;
151:   Vec            localU;

154:   TSGetDM(ts,&da);
155:   DMGetLocalVector(da,&localU);
156:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

158:   hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx);
159:   hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy);

161:   /*
162:      Scatter ghost points to local vector,using the 2-step process
163:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
164:      By placing code between these two statements, computations can be
165:      done while messages are in transition.
166:   */
167:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
168:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

170:   /*
171:      Get pointers to vector data
172:   */
173:   DMDAVecGetArrayRead(da,localU,&u);
174:   DMDAVecGetArray(da,F,&f);

176:   /*
177:      Get local grid boundaries
178:   */
179:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

181:   /*
182:      Compute function over the locally owned part of the grid
183:   */
184:   for (j=ys; j<ys+ym; j++) {
185:     for (i=xs; i<xs+xm; i++) {
186:       uc        = u[j][i].u;
187:       uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
188:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
189:       vc        = u[j][i].v;
190:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
191:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
192:       f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);
193:       f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc;
194:     }
195:   }
196:   PetscLogFlops(16*xm*ym);

198:   /*
199:      Restore vectors
200:   */
201:   DMDAVecRestoreArrayRead(da,localU,&u);
202:   DMDAVecRestoreArray(da,F,&f);
203:   DMRestoreLocalVector(da,&localU);
204:   return(0);
205: }

207: /* ------------------------------------------------------------------- */
208: PetscErrorCode InitialConditions(DM da,Vec U)
209: {
211:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
212:   Field          **u;
213:   PetscReal      hx,hy,x,y;

216:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

218:   hx = 2.5/(PetscReal)(Mx);
219:   hy = 2.5/(PetscReal)(My);

221:   /*
222:      Get pointers to vector data
223:   */
224:   DMDAVecGetArray(da,U,&u);

226:   /*
227:      Get local grid boundaries
228:   */
229:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

231:   /*
232:      Compute function over the locally owned part of the grid
233:   */
234:   for (j=ys; j<ys+ym; j++) {
235:     y = j*hy;
236:     for (i=xs; i<xs+xm; i++) {
237:       x = i*hx;
238:       if ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0);
239:       else u[j][i].v = 0.0;

241:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
242:     }
243:   }

245:   /*
246:      Restore vectors
247:   */
248:   DMDAVecRestoreArray(da,U,&u);
249:   return(0);
250: }

252: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
253: {
254:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
255:   DM             da;
257:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
258:   PetscReal      hx,hy,sx,sy;
259:   PetscScalar    uc,vc;
260:   Field          **u;
261:   Vec            localU;
262:   MatStencil     stencil[6],rowstencil;
263:   PetscScalar    entries[6];

266:   TSGetDM(ts,&da);
267:   DMGetLocalVector(da,&localU);
268:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

270:   hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx);
271:   hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy);

273:   /*
274:      Scatter ghost points to local vector,using the 2-step process
275:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
276:      By placing code between these two statements, computations can be
277:      done while messages are in transition.
278:   */
279:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
280:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

282:   /*
283:      Get pointers to vector data
284:   */
285:   DMDAVecGetArrayRead(da,localU,&u);

287:   /*
288:      Get local grid boundaries
289:   */
290:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

292:   stencil[0].k = 0;
293:   stencil[1].k = 0;
294:   stencil[2].k = 0;
295:   stencil[3].k = 0;
296:   stencil[4].k = 0;
297:   stencil[5].k = 0;
298:   rowstencil.k = 0;
299:   /*
300:      Compute function over the locally owned part of the grid
301:   */
302:   for (j=ys; j<ys+ym; j++) {

304:     stencil[0].j = j-1;
305:     stencil[1].j = j+1;
306:     stencil[2].j = j;
307:     stencil[3].j = j;
308:     stencil[4].j = j;
309:     stencil[5].j = j;
310:     rowstencil.k = 0; rowstencil.j = j;
311:     for (i=xs; i<xs+xm; i++) {
312:       uc = u[j][i].u;
313:       vc = u[j][i].v;

315:       /*      uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
316:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;

318:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
319:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
320:        f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/

322:       stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy;
323:       stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy;
324:       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx;
325:       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx;
326:       stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma;
327:       stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc;
328:       rowstencil.i = i; rowstencil.c = 0;

330:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);

332:       stencil[0].c = 1; entries[0] = appctx->D2*sy;
333:       stencil[1].c = 1; entries[1] = appctx->D2*sy;
334:       stencil[2].c = 1; entries[2] = appctx->D2*sx;
335:       stencil[3].c = 1; entries[3] = appctx->D2*sx;
336:       stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa;
337:       stencil[5].c = 0; entries[5] = vc*vc;
338:       rowstencil.c = 1;

340:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
341:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
342:     }
343:   }

345:   /*
346:      Restore vectors
347:   */
348:   PetscLogFlops(19*xm*ym);
349:   DMDAVecRestoreArrayRead(da,localU,&u);
350:   DMRestoreLocalVector(da,&localU);
351:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
352:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
353:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
354:   return(0);
355: }


358: /*TEST

360:    test:
361:       args: -ts_view  -ts_monitor -ts_max_time 500
362:       requires: double
363:       timeoutfactor: 3

365:    test:
366:       suffix: 2
367:       args: -ts_view  -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution
368:       requires: x double
369:       output_file: output/ex5_1.out
370:       timeoutfactor: 3

372: TEST*/