Actual source code: ex5.c
1: static char help[] = "Demonstrates Pattern Formation with Reaction-Diffusion Equations.\n";
3: /*F
4: This example is taken from the book, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
5: W. Hundsdorf and J.G. Verwer, Page 21, Pattern Formation with Reaction-Diffusion Equations
6: \begin{eqnarray*}
7: u_t = D_1 (u_{xx} + u_{yy}) - u*v^2 + \gamma(1 -u) \\
8: v_t = D_2 (v_{xx} + v_{yy}) + u*v^2 - (\gamma + \kappa)v
9: \end{eqnarray*}
10: Unlike in the book this uses periodic boundary conditions instead of Neumann
11: (since they are easier for finite differences).
12: F*/
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 numerical (ODE) integrators. Note that this
31: file automatically includes:
32: petscsys.h - base PETSc routines
33: petscis.h - index sets
34: petscvec.h - vectors
35: petscmat.h - matrices
36: petscpc.h - preconditioners
37: petscksp.h - Krylov subspace methods
38: petscsnes.h - nonlinear solvers
39: */
40: #include "reaction_diffusion.h"
41: #include <petscdm.h>
42: #include <petscdmda.h>
44: /* ------------------------------------------------------------------- */
45: PetscErrorCode InitialConditions(DM da, Vec U)
46: {
47: PetscInt i, j, xs, ys, xm, ym, Mx, My;
48: Field **u;
49: PetscReal hx, hy, x, y;
51: PetscFunctionBegin;
52: PetscCall(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));
54: hx = 2.5 / (PetscReal)Mx;
55: hy = 2.5 / (PetscReal)My;
57: /*
58: Get pointers to actual vector data
59: */
60: PetscCall(DMDAVecGetArray(da, U, &u));
62: /*
63: Get local grid boundaries
64: */
65: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
67: /*
68: Compute function over the locally owned part of the grid
69: */
70: for (j = ys; j < ys + ym; j++) {
71: y = j * hy;
72: for (i = xs; i < xs + xm; i++) {
73: x = i * hx;
74: if (PetscApproximateGTE(x, 1.0) && PetscApproximateLTE(x, 1.5) && PetscApproximateGTE(y, 1.0) && PetscApproximateLTE(y, 1.5))
75: u[j][i].v = PetscPowReal(PetscSinReal(4.0 * PETSC_PI * x), 2.0) * PetscPowReal(PetscSinReal(4.0 * PETSC_PI * y), 2.0) / 4.0;
76: else u[j][i].v = 0.0;
78: u[j][i].u = 1.0 - 2.0 * u[j][i].v;
79: }
80: }
82: /*
83: Restore access to vector
84: */
85: PetscCall(DMDAVecRestoreArray(da, U, &u));
86: PetscFunctionReturn(PETSC_SUCCESS);
87: }
89: int main(int argc, char **argv)
90: {
91: TS ts; /* ODE integrator */
92: Vec x; /* solution */
93: DM da;
94: AppCtx appctx;
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Initialize program
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscFunctionBeginUser;
100: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
101: appctx.D1 = 8.0e-5;
102: appctx.D2 = 4.0e-5;
103: appctx.gamma = .024;
104: appctx.kappa = .06;
105: appctx.aijpc = PETSC_FALSE;
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create distributed array (DMDA) to manage parallel grid and vectors
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 65, 65, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, NULL, &da));
111: PetscCall(DMSetFromOptions(da));
112: PetscCall(DMSetUp(da));
113: PetscCall(DMDASetFieldName(da, 0, "u"));
114: PetscCall(DMDASetFieldName(da, 1, "v"));
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Create global vector from DMDA; this will be used to store the solution
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: PetscCall(DMCreateGlobalVector(da, &x));
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create timestepping solver context
123: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
125: PetscCall(TSSetType(ts, TSARKIMEX));
126: PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE));
127: PetscCall(TSSetDM(ts, da));
128: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
129: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &appctx));
130: PetscCall(TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, &appctx));
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Set initial conditions
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135: PetscCall(InitialConditions(da, x));
136: PetscCall(TSSetSolution(ts, x));
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Set solver options
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: PetscCall(TSSetMaxTime(ts, 2000.0));
142: PetscCall(TSSetTimeStep(ts, .0001));
143: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
144: PetscCall(TSSetFromOptions(ts));
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Solve ODE system
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: PetscCall(TSSolve(ts, x));
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Free work space.
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: PetscCall(VecDestroy(&x));
155: PetscCall(TSDestroy(&ts));
156: PetscCall(DMDestroy(&da));
158: PetscCall(PetscFinalize());
159: return 0;
160: }
162: /*TEST
164: build:
165: depends: reaction_diffusion.c
167: test:
168: args: -ts_view -ts_monitor -ts_max_time 500
169: requires: double
170: timeoutfactor: 3
172: test:
173: suffix: 2
174: args: -ts_view -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution
175: requires: x double
176: output_file: output/ex5_1.out
177: timeoutfactor: 3
179: TEST*/