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 <petscdm.h> 39 #include <petscdmda.h> 40 #include <petscts.h> 41 42 /* Simple C struct that allows us to access the two velocity (x and y directions) values easily in the code */ 43 typedef struct { 44 PetscScalar u,v; 45 } Field; 46 47 /* Data structure to store the model parameters */ 48 typedef struct { 49 PetscReal D1,D2,gamma,kappa; 50 } AppCtx; 51 52 /* 53 User-defined routines 54 */ 55 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*),InitialConditions(DM,Vec); 56 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 57 58 int main(int argc,char **argv) 59 { 60 TS ts; /* ODE integrator */ 61 Vec x; /* solution */ 62 PetscErrorCode ierr; 63 DM da; 64 AppCtx appctx; 65 66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 67 Initialize program 68 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 69 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 70 PetscFunctionBeginUser; 71 appctx.D1 = 8.0e-5; 72 appctx.D2 = 4.0e-5; 73 appctx.gamma = .024; 74 appctx.kappa = .06; 75 76 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 77 Create distributed array (DMDA) to manage parallel grid and vectors 78 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 79 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,65,65,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr); 80 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 81 ierr = DMSetUp(da);CHKERRQ(ierr); 82 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr); 83 ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr); 84 85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 86 Create global vector from DMDA; this will be used to store the solution 87 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 88 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr); 89 90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 91 Create timestepping solver context 92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 93 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 94 ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr); 95 ierr = TSARKIMEXSetFullyImplicit(ts,PETSC_TRUE);CHKERRQ(ierr); 96 ierr = TSSetDM(ts,da);CHKERRQ(ierr); 97 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 98 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);CHKERRQ(ierr); 99 ierr = TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr); 100 101 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 102 Set initial conditions 103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 104 ierr = InitialConditions(da,x);CHKERRQ(ierr); 105 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 106 107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 108 Set solver options 109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 110 ierr = TSSetMaxTime(ts,2000.0);CHKERRQ(ierr); 111 ierr = TSSetTimeStep(ts,.0001);CHKERRQ(ierr); 112 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 113 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 114 115 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 116 Solve ODE system 117 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 118 ierr = TSSolve(ts,x);CHKERRQ(ierr); 119 120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 121 Free work space. 122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 123 ierr = VecDestroy(&x);CHKERRQ(ierr); 124 ierr = TSDestroy(&ts);CHKERRQ(ierr); 125 ierr = DMDestroy(&da);CHKERRQ(ierr); 126 127 ierr = PetscFinalize(); 128 return ierr; 129 } 130 /* ------------------------------------------------------------------- */ 131 /* 132 RHSFunction - Evaluates nonlinear function, that defines the right 133 hand side of the ODE 134 135 Input Parameters: 136 . ts - the TS context 137 . time - current time 138 . U - input vector 139 . ptr - optional user-defined context, as set by TSSetRHSFunction() 140 141 Output Parameter: 142 . F - function vector 143 */ 144 PetscErrorCode RHSFunction(TS ts,PetscReal time,Vec U,Vec F,void *ptr) 145 { 146 AppCtx *appctx = (AppCtx*)ptr; 147 DM da; 148 PetscErrorCode ierr; 149 PetscInt i,j,Mx,My,xs,ys,xm,ym; 150 PetscReal hx,hy,sx,sy; 151 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 152 Field **u,**f; 153 Vec localU; 154 155 PetscFunctionBegin; 156 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 157 /* Get local (ghosted) work vector */ 158 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 159 /* Get information about mesh needed for discretization */ 160 ierr = 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);CHKERRQ(ierr); 161 162 hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx); 163 hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy); 164 165 /* 166 Scatter ghost points to local vector, using the 2-step process 167 */ 168 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 169 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 170 171 /* 172 Get pointers to actual vector data 173 */ 174 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 175 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 176 177 /* 178 Get local grid boundaries; this is the region that this process owns and must operate on 179 */ 180 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 181 182 /* 183 Compute function over the locally owned part of the grid with standard finite differences 184 */ 185 for (j=ys; j<ys+ym; j++) { 186 for (i=xs; i<xs+xm; i++) { 187 uc = u[j][i].u; 188 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 189 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 190 vc = u[j][i].v; 191 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 192 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 193 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); 194 f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; 195 } 196 } 197 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 198 199 /* 200 Restore access to vectors and return no longer needed work vector 201 */ 202 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 203 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 204 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 205 PetscFunctionReturn(0); 206 } 207 208 /* ------------------------------------------------------------------- */ 209 PetscErrorCode InitialConditions(DM da,Vec U) 210 { 211 PetscErrorCode ierr; 212 PetscInt i,j,xs,ys,xm,ym,Mx,My; 213 Field **u; 214 PetscReal hx,hy,x,y; 215 216 PetscFunctionBegin; 217 ierr = 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);CHKERRQ(ierr); 218 219 hx = 2.5/(PetscReal)(Mx); 220 hy = 2.5/(PetscReal)(My); 221 222 /* 223 Get pointers to actual vector data 224 */ 225 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 226 227 /* 228 Get local grid boundaries 229 */ 230 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 231 232 /* 233 Compute function over the locally owned part of the grid 234 */ 235 for (j=ys; j<ys+ym; j++) { 236 y = j*hy; 237 for (i=xs; i<xs+xm; i++) { 238 x = i*hx; 239 if (PetscGTE(x,1.0) && PetscLTE(x,1.5) && PetscGTE(y,1.0) && PetscLTE(y,1.5)) u[j][i].v = PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0)/4.0; 240 else u[j][i].v = 0.0; 241 242 u[j][i].u = 1.0 - 2.0*u[j][i].v; 243 } 244 } 245 246 /* 247 Restore access to vector 248 */ 249 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 /* 254 RHSJacobian - Evaluates the Jacobian of the right hand side 255 function of the ODE. 256 257 Input Parameters: 258 . ts - the TS context 259 . time - current time 260 . U - input vector 261 . ptr - optional user-defined context, as set by TSSetRHSJacobian() 262 263 Output Parameter: 264 . A - the Jacobian 265 . BB - optional additional matrix where an approximation to the Jacobian 266 may be stored from which the preconditioner is constructed 267 */ 268 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx) 269 { 270 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 271 DM da; 272 PetscErrorCode ierr; 273 PetscInt i,j,Mx,My,xs,ys,xm,ym; 274 PetscReal hx,hy,sx,sy; 275 PetscScalar uc,vc; 276 Field **u; 277 Vec localU; 278 MatStencil stencil[6],rowstencil; 279 PetscScalar entries[6]; 280 281 PetscFunctionBegin; 282 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 283 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 284 ierr = 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);CHKERRQ(ierr); 285 286 hx = 2.50/(PetscReal)(Mx); sx = 1.0/(hx*hx); 287 hy = 2.50/(PetscReal)(My); sy = 1.0/(hy*hy); 288 289 /* 290 Scatter ghost points to local vector,using the 2-step process 291 */ 292 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 293 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 294 295 /* 296 Get pointers to vector data 297 */ 298 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 299 300 /* 301 Get local grid boundaries 302 */ 303 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 304 305 stencil[0].k = 0; 306 stencil[1].k = 0; 307 stencil[2].k = 0; 308 stencil[3].k = 0; 309 stencil[4].k = 0; 310 stencil[5].k = 0; 311 rowstencil.k = 0; 312 /* 313 Compute function over the locally owned part of the grid 314 */ 315 for (j=ys; j<ys+ym; j++) { 316 317 stencil[0].j = j-1; 318 stencil[1].j = j+1; 319 stencil[2].j = j; 320 stencil[3].j = j; 321 stencil[4].j = j; 322 stencil[5].j = j; 323 rowstencil.k = 0; rowstencil.j = j; 324 for (i=xs; i<xs+xm; i++) { 325 uc = u[j][i].u; 326 vc = u[j][i].v; 327 328 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 329 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 330 331 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 332 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 333 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/ 334 335 stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy; 336 stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy; 337 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx; 338 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx; 339 stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma; 340 stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc; 341 rowstencil.i = i; rowstencil.c = 0; 342 343 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 344 345 stencil[0].c = 1; entries[0] = appctx->D2*sy; 346 stencil[1].c = 1; entries[1] = appctx->D2*sy; 347 stencil[2].c = 1; entries[2] = appctx->D2*sx; 348 stencil[3].c = 1; entries[3] = appctx->D2*sx; 349 stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa; 350 stencil[5].c = 0; entries[5] = vc*vc; 351 rowstencil.c = 1; 352 353 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 354 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 355 } 356 } 357 358 /* 359 Restore vectors 360 */ 361 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 362 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 363 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 364 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 365 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 366 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 367 PetscFunctionReturn(0); 368 } 369 370 /*TEST 371 372 test: 373 args: -ts_view -ts_monitor -ts_max_time 500 374 requires: double 375 timeoutfactor: 3 376 377 test: 378 suffix: 2 379 args: -ts_view -ts_monitor -ts_max_time 500 -ts_monitor_draw_solution 380 requires: x double 381 output_file: output/ex5_1.out 382 timeoutfactor: 3 383 384 TEST*/ 385