1 static char help[] = "Demonstrates Pattern Formation with Reaction-Diffusion Equations.\n"; 2 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*/ 13 14 /* 15 Helpful runtime monitor options: 16 -ts_monitor_draw_solution 17 -draw_save -draw_save_movie 18 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) 21 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 . 24 25 */ 26 27 /* 28 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 petscvec.h - vectors 33 petscmat.h - matrices petscis.h - index sets 34 petscksp.h - Krylov subspace methods petscpc.h - preconditioners 35 petscviewer.h - viewers petscsnes.h - nonlinear solvers 36 */ 37 #include "reaction_diffusion.h" 38 #include <petscdm.h> 39 #include <petscdmda.h> 40 41 /* ------------------------------------------------------------------- */ 42 PetscErrorCode InitialConditions(DM da, Vec U) 43 { 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(PETSC_SUCCESS); 84 } 85 86 int main(int argc, char **argv) 87 { 88 TS ts; /* ODE integrator */ 89 Vec x; /* solution */ 90 DM da; 91 AppCtx appctx; 92 93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 94 Initialize program 95 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 96 PetscFunctionBeginUser; 97 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 98 PetscFunctionBeginUser; 99 appctx.D1 = 8.0e-5; 100 appctx.D2 = 4.0e-5; 101 appctx.gamma = .024; 102 appctx.kappa = .06; 103 appctx.aijpc = PETSC_FALSE; 104 105 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 106 Create distributed array (DMDA) to manage parallel grid and vectors 107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 108 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)); 109 PetscCall(DMSetFromOptions(da)); 110 PetscCall(DMSetUp(da)); 111 PetscCall(DMDASetFieldName(da, 0, "u")); 112 PetscCall(DMDASetFieldName(da, 1, "v")); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Create global vector from DMDA; this will be used to store the solution 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 PetscCall(DMCreateGlobalVector(da, &x)); 118 119 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 120 Create timestepping solver context 121 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 122 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 123 PetscCall(TSSetType(ts, TSARKIMEX)); 124 PetscCall(TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE)); 125 PetscCall(TSSetDM(ts, da)); 126 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 127 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &appctx)); 128 PetscCall(TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, &appctx)); 129 130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 131 Set initial conditions 132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 133 PetscCall(InitialConditions(da, x)); 134 PetscCall(TSSetSolution(ts, x)); 135 136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 137 Set solver options 138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 139 PetscCall(TSSetMaxTime(ts, 2000.0)); 140 PetscCall(TSSetTimeStep(ts, .0001)); 141 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 142 PetscCall(TSSetFromOptions(ts)); 143 144 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 145 Solve ODE system 146 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 147 PetscCall(TSSolve(ts, x)); 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Free work space. 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 PetscCall(VecDestroy(&x)); 153 PetscCall(TSDestroy(&ts)); 154 PetscCall(DMDestroy(&da)); 155 156 PetscCall(PetscFinalize()); 157 return 0; 158 } 159 160 /*TEST 161 162 build: 163 depends: reaction_diffusion.c 164 165 test: 166 args: -ts_view -ts_monitor -ts_max_time 500 167 requires: double 168 timeoutfactor: 3 169 170 test: 171 suffix: 2 172 args: -ts_view -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution 173 requires: x double 174 output_file: output/ex5_1.out 175 timeoutfactor: 3 176 177 TEST*/ 178