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 PetscInt i, j, xs, ys, xm, ym, Mx, My; 45 Field **u; 46 PetscReal hx, hy, x, y; 47 48 PetscFunctionBegin; 49 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)); 50 51 hx = 2.5 / (PetscReal)(Mx); 52 hy = 2.5 / (PetscReal)(My); 53 54 /* 55 Get pointers to actual vector data 56 */ 57 PetscCall(DMDAVecGetArray(da, U, &u)); 58 59 /* 60 Get local grid boundaries 61 */ 62 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 63 64 /* 65 Compute function over the locally owned part of the grid 66 */ 67 for (j = ys; j < ys + ym; j++) { 68 y = j * hy; 69 for (i = xs; i < xs + xm; i++) { 70 x = i * hx; 71 if (PetscApproximateGTE(x, 1.0) && PetscApproximateLTE(x, 1.5) && PetscApproximateGTE(y, 1.0) && PetscApproximateLTE(y, 1.5)) 72 u[j][i].v = PetscPowReal(PetscSinReal(4.0 * PETSC_PI * x), 2.0) * PetscPowReal(PetscSinReal(4.0 * PETSC_PI * y), 2.0) / 4.0; 73 else u[j][i].v = 0.0; 74 75 u[j][i].u = 1.0 - 2.0 * u[j][i].v; 76 } 77 } 78 79 /* 80 Restore access to vector 81 */ 82 PetscCall(DMDAVecRestoreArray(da, U, &u)); 83 PetscFunctionReturn(0); 84 } 85 86 int main(int argc, char **argv) { 87 TS ts; /* ODE integrator */ 88 Vec x; /* solution */ 89 DM da; 90 AppCtx appctx; 91 92 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 93 Initialize program 94 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 95 PetscFunctionBeginUser; 96 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 97 PetscFunctionBeginUser; 98 appctx.D1 = 8.0e-5; 99 appctx.D2 = 4.0e-5; 100 appctx.gamma = .024; 101 appctx.kappa = .06; 102 appctx.aijpc = PETSC_FALSE; 103 104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 Create distributed array (DMDA) to manage parallel grid and vectors 106 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 107 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)); 108 PetscCall(DMSetFromOptions(da)); 109 PetscCall(DMSetUp(da)); 110 PetscCall(DMDASetFieldName(da, 0, "u")); 111 PetscCall(DMDASetFieldName(da, 1, "v")); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Create global vector from DMDA; this will be used to store the solution 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 PetscCall(DMCreateGlobalVector(da, &x)); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Create timestepping solver context 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 122 PetscCall(TSSetType(ts, TSARKIMEX)); 123 PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE)); 124 PetscCall(TSSetDM(ts, da)); 125 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 126 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &appctx)); 127 PetscCall(TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, &appctx)); 128 129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 Set initial conditions 131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 132 PetscCall(InitialConditions(da, x)); 133 PetscCall(TSSetSolution(ts, x)); 134 135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 Set solver options 137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 138 PetscCall(TSSetMaxTime(ts, 2000.0)); 139 PetscCall(TSSetTimeStep(ts, .0001)); 140 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 141 PetscCall(TSSetFromOptions(ts)); 142 143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 144 Solve ODE system 145 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 146 PetscCall(TSSolve(ts, x)); 147 148 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 149 Free work space. 150 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 151 PetscCall(VecDestroy(&x)); 152 PetscCall(TSDestroy(&ts)); 153 PetscCall(DMDestroy(&da)); 154 155 PetscCall(PetscFinalize()); 156 return 0; 157 } 158 159 /*TEST 160 161 build: 162 depends: reaction_diffusion.c 163 164 test: 165 args: -ts_view -ts_monitor -ts_max_time 500 166 requires: double 167 timeoutfactor: 3 168 169 test: 170 suffix: 2 171 args: -ts_view -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution 172 requires: x double 173 output_file: output/ex5_1.out 174 timeoutfactor: 3 175 176 TEST*/ 177