Actual source code: ex5.c
petsc-3.10.3 2018-12-18
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*/