Actual source code: ex5.c

petsc-3.12.2 2019-11-22
Report Typos and Errors

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

This example is taken from the book, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
W. Hundsdorf and J.G. Verwer, Page 21, Pattern Formation with Reaction-Diffusion Equations
\begin{eqnarray*}
u_t = D_1 (u_{xx} + u_{yy}) - u*v^2 + \gamma(1 -u) \\
v_t = D_2 (v_{xx} + v_{yy}) + u*v^2 - (\gamma + \kappa)v
\end{eqnarray*}
Unlike in the book this uses periodic boundary conditions instead of Neumann
(since they are easier for finite differences).
 15: /*
 16:       Helpful runtime monitor options:
 17:            -ts_monitor_draw_solution
 18:            -draw_save -draw_save_movie

 20:       Helpful runtime linear solver options:
 21:            -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)

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

 26: */

 28: /*

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

 42: /*  Simple C struct that allows us to access the two velocity (x and y directions) values easily in the code */
 43: typedef struct {
 44:   PetscScalar u,v;
 45: } Field;

 47: /*  Data structure to store the model parameters */
 48: typedef struct {
 49:   PetscReal D1,D2,gamma,kappa;
 50: } AppCtx;

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

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

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

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

 85:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 86:      Create global vector from DMDA; this will be used to store the solution
 87:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 88:   DMCreateGlobalVector(da,&x);

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

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

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

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

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121:      Free work space.  
122:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123:   VecDestroy(&x);
124:   TSDestroy(&ts);
125:   DMDestroy(&da);

127:   PetscFinalize();
128:   return ierr;
129: }
130: /* ------------------------------------------------------------------- */
131: /*
132:    RHSFunction - Evaluates nonlinear function, that defines the right 
133:      hand side of the ODE

135:    Input Parameters:
136: .  ts - the TS context
137: .  time - current time
138: .  U - input vector
139: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

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

156:   TSGetDM(ts,&da);
157:   /* Get local (ghosted) work vector */
158:   DMGetLocalVector(da,&localU);
159:   /* Get information about mesh needed for discretization */
160:   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);

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

165:   /*
166:      Scatter ghost points to local vector, using the 2-step process
167:   */
168:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
169:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

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

177:   /*
178:      Get local grid boundaries; this is the region that this process owns and must operate on
179:   */
180:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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

199:   /*
200:      Restore access to vectors and return no longer needed work vector
201:   */
202:   DMDAVecRestoreArrayRead(da,localU,&u);
203:   DMDAVecRestoreArray(da,F,&f);
204:   DMRestoreLocalVector(da,&localU);
205:   return(0);
206: }

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

217:   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);

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

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

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

232:   /*
233:      Compute function over the locally owned part of the grid
234:   */
235:   for (j=ys; j<ys+ym; j++) {
236:     y = j*hy;
237:     for (i=xs; i<xs+xm; i++) {
238:       x = i*hx;
239:       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);
240:       else u[j][i].v = 0.0;

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

246:   /*
247:      Restore access to vector
248:   */
249:   DMDAVecRestoreArray(da,U,&u);
250:   return(0);
251: }

253: /*
254:    RHSJacobian - Evaluates the Jacobian of the right hand side
255:      function of the ODE.

257:    Input Parameters:
258: .  ts - the TS context
259: .  time - current time
260: .  U - input vector
261: .  ptr - optional user-defined context, as set by TSSetRHSJacobian()

263:    Output Parameter:
264: .  A - the Jacobian
265: .  BB - optional additional matrix where an approximation to the Jacobian
266:         may be stored from which the preconditioner is constructed
267:  */
268: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
269: {
270:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
271:   DM             da;
273:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
274:   PetscReal      hx,hy,sx,sy;
275:   PetscScalar    uc,vc;
276:   Field          **u;
277:   Vec            localU;
278:   MatStencil     stencil[6],rowstencil;
279:   PetscScalar    entries[6];

282:   TSGetDM(ts,&da);
283:   DMGetLocalVector(da,&localU);
284:   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);

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

289:   /*
290:      Scatter ghost points to local vector,using the 2-step process
291:   */
292:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
293:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

295:   /*
296:      Get pointers to vector data
297:   */
298:   DMDAVecGetArrayRead(da,localU,&u);

300:   /*
301:      Get local grid boundaries
302:   */
303:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

305:   stencil[0].k = 0;
306:   stencil[1].k = 0;
307:   stencil[2].k = 0;
308:   stencil[3].k = 0;
309:   stencil[4].k = 0;
310:   stencil[5].k = 0;
311:   rowstencil.k = 0;
312:   /*
313:      Compute function over the locally owned part of the grid
314:   */
315:   for (j=ys; j<ys+ym; j++) {

317:     stencil[0].j = j-1;
318:     stencil[1].j = j+1;
319:     stencil[2].j = j;
320:     stencil[3].j = j;
321:     stencil[4].j = j;
322:     stencil[5].j = j;
323:     rowstencil.k = 0; rowstencil.j = j;
324:     for (i=xs; i<xs+xm; i++) {
325:       uc = u[j][i].u;
326:       vc = u[j][i].v;

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

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

335:       stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy;
336:       stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy;
337:       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx;
338:       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx;
339:       stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma;
340:       stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc;
341:       rowstencil.i = i; rowstencil.c = 0;

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

345:       stencil[0].c = 1; entries[0] = appctx->D2*sy;
346:       stencil[1].c = 1; entries[1] = appctx->D2*sy;
347:       stencil[2].c = 1; entries[2] = appctx->D2*sx;
348:       stencil[3].c = 1; entries[3] = appctx->D2*sx;
349:       stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa;
350:       stencil[5].c = 0; entries[5] = vc*vc;
351:       rowstencil.c = 1;

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

358:   /*
359:      Restore vectors
360:   */
361:   PetscLogFlops(19*xm*ym);
362:   DMDAVecRestoreArrayRead(da,localU,&u);
363:   DMRestoreLocalVector(da,&localU);
364:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
365:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
366:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
367:   return(0);
368: }


371: /*TEST

373:    test:
374:       args: -ts_view  -ts_monitor -ts_max_time 500
375:       requires: double
376:       timeoutfactor: 3

378:    test:
379:       suffix: 2
380:       args: -ts_view  -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution
381:       requires: x double
382:       output_file: output/ex5_1.out
383:       timeoutfactor: 3

385: TEST*/