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