1 #include "reaction_diffusion.h" 2 #include <petscdm.h> 3 #include <petscdmda.h> 4 5 /*F 6 This example is taken from the book, Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by 7 W. Hundsdorf and J.G. Verwer, Page 21, Pattern Formation with Reaction-Diffusion Equations 8 \begin{eqnarray*} 9 u_t = D_1 (u_{xx} + u_{yy}) - u*v^2 + \gamma(1 -u) \\ 10 v_t = D_2 (v_{xx} + v_{yy}) + u*v^2 - (\gamma + \kappa)v 11 \end{eqnarray*} 12 Unlike in the book this uses periodic boundary conditions instead of Neumann 13 (since they are easier for finite differences). 14 F*/ 15 16 /* 17 RHSFunction - Evaluates nonlinear function, F(x). 18 19 Input Parameters: 20 . ts - the TS context 21 . X - input vector 22 . ptr - optional user-defined context, as set by TSSetRHSFunction() 23 24 Output Parameter: 25 . F - function vector 26 */ 27 PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr) 28 { 29 AppCtx *appctx = (AppCtx*)ptr; 30 DM da; 31 PetscErrorCode ierr; 32 PetscInt i,j,Mx,My,xs,ys,xm,ym; 33 PetscReal hx,hy,sx,sy; 34 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 35 Field **u,**f; 36 Vec localU; 37 38 PetscFunctionBegin; 39 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 40 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 41 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); 42 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 43 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 44 45 /* 46 Scatter ghost points to local vector,using the 2-step process 47 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 48 By placing code between these two statements, computations can be 49 done while messages are in transition. 50 */ 51 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 52 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 53 54 /* 55 Get pointers to vector data 56 */ 57 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 58 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 59 60 /* 61 Get local grid boundaries 62 */ 63 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 64 65 /* 66 Compute function over the locally owned part of the grid 67 */ 68 for (j=ys; j<ys+ym; j++) { 69 for (i=xs; i<xs+xm; i++) { 70 uc = u[j][i].u; 71 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 72 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 73 vc = u[j][i].v; 74 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 75 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 76 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); 77 f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; 78 } 79 } 80 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 81 82 /* 83 Restore vectors 84 */ 85 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 86 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 87 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 88 PetscFunctionReturn(0); 89 } 90 91 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx) 92 { 93 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 94 DM da; 95 PetscErrorCode ierr; 96 PetscInt i,j,Mx,My,xs,ys,xm,ym; 97 PetscReal hx,hy,sx,sy; 98 PetscScalar uc,vc; 99 Field **u; 100 Vec localU; 101 MatStencil stencil[6],rowstencil; 102 PetscScalar entries[6]; 103 104 PetscFunctionBegin; 105 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 106 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 107 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); 108 109 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 110 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 111 112 /* 113 Scatter ghost points to local vector,using the 2-step process 114 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 115 By placing code between these two statements, computations can be 116 done while messages are in transition. 117 */ 118 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 119 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 120 121 /* 122 Get pointers to vector data 123 */ 124 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 125 126 /* 127 Get local grid boundaries 128 */ 129 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 130 131 stencil[0].k = 0; 132 stencil[1].k = 0; 133 stencil[2].k = 0; 134 stencil[3].k = 0; 135 stencil[4].k = 0; 136 stencil[5].k = 0; 137 rowstencil.k = 0; 138 /* 139 Compute function over the locally owned part of the grid 140 */ 141 for (j=ys; j<ys+ym; j++) { 142 stencil[0].j = j-1; 143 stencil[1].j = j+1; 144 stencil[2].j = j; 145 stencil[3].j = j; 146 stencil[4].j = j; 147 stencil[5].j = j; 148 rowstencil.k = 0; rowstencil.j = j; 149 for (i=xs; i<xs+xm; i++) { 150 uc = u[j][i].u; 151 vc = u[j][i].v; 152 153 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 154 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 155 156 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 157 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 158 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/ 159 160 stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy; 161 stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy; 162 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx; 163 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx; 164 stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma; 165 stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc; 166 rowstencil.i = i; rowstencil.c = 0; 167 168 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 169 if (appctx->aijpc) { 170 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 171 } 172 stencil[0].c = 1; entries[0] = appctx->D2*sy; 173 stencil[1].c = 1; entries[1] = appctx->D2*sy; 174 stencil[2].c = 1; entries[2] = appctx->D2*sx; 175 stencil[3].c = 1; entries[3] = appctx->D2*sx; 176 stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa; 177 stencil[5].c = 0; entries[5] = vc*vc; 178 rowstencil.c = 1; 179 180 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 181 if (appctx->aijpc) { 182 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 183 } 184 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 185 } 186 } 187 188 /* 189 Restore vectors 190 */ 191 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 192 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 193 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 194 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 195 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 196 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 197 if (appctx->aijpc) { 198 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 199 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 200 ierr = MatSetOption(BB,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 201 } 202 PetscFunctionReturn(0); 203 } 204 205 /* 206 IFunction - Evaluates implicit nonlinear function, xdot - F(x). 207 208 Input Parameters: 209 . ts - the TS context 210 . U - input vector 211 . Udot - input vector 212 . ptr - optional user-defined context, as set by TSSetRHSFunction() 213 214 Output Parameter: 215 . F - function vector 216 */ 217 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr) 218 { 219 AppCtx *appctx = (AppCtx*)ptr; 220 DM da; 221 PetscErrorCode ierr; 222 PetscInt i,j,Mx,My,xs,ys,xm,ym; 223 PetscReal hx,hy,sx,sy; 224 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 225 Field **u,**f,**udot; 226 Vec localU; 227 228 PetscFunctionBegin; 229 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 230 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 231 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); 232 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 233 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 234 235 /* 236 Scatter ghost points to local vector,using the 2-step process 237 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 238 By placing code between these two statements, computations can be 239 done while messages are in transition. 240 */ 241 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 242 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 243 244 /* 245 Get pointers to vector data 246 */ 247 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 248 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 249 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr); 250 251 /* 252 Get local grid boundaries 253 */ 254 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 255 256 /* 257 Compute function over the locally owned part of the grid 258 */ 259 for (j=ys; j<ys+ym; j++) { 260 for (i=xs; i<xs+xm; i++) { 261 uc = u[j][i].u; 262 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 263 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 264 vc = u[j][i].v; 265 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 266 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 267 f[j][i].u = udot[j][i].u - ( appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc)); 268 f[j][i].v = udot[j][i].v - ( appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc); 269 } 270 } 271 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 272 273 /* 274 Restore vectors 275 */ 276 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 277 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 278 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr); 279 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void *ctx) 284 { 285 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 286 DM da; 287 PetscErrorCode ierr; 288 PetscInt i,j,Mx,My,xs,ys,xm,ym; 289 PetscReal hx,hy,sx,sy; 290 PetscScalar uc,vc; 291 Field **u; 292 Vec localU; 293 MatStencil stencil[6],rowstencil; 294 PetscScalar entries[6]; 295 296 PetscFunctionBegin; 297 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 298 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 299 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); 300 301 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 302 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 303 304 /* 305 Scatter ghost points to local vector,using the 2-step process 306 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 307 By placing code between these two statements, computations can be 308 done while messages are in transition. 309 */ 310 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 311 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 312 313 /* 314 Get pointers to vector data 315 */ 316 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 317 318 /* 319 Get local grid boundaries 320 */ 321 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 322 323 stencil[0].k = 0; 324 stencil[1].k = 0; 325 stencil[2].k = 0; 326 stencil[3].k = 0; 327 stencil[4].k = 0; 328 stencil[5].k = 0; 329 rowstencil.k = 0; 330 /* 331 Compute function over the locally owned part of the grid 332 */ 333 for (j=ys; j<ys+ym; j++) { 334 335 stencil[0].j = j-1; 336 stencil[1].j = j+1; 337 stencil[2].j = j; 338 stencil[3].j = j; 339 stencil[4].j = j; 340 stencil[5].j = j; 341 rowstencil.k = 0; rowstencil.j = j; 342 for (i=xs; i<xs+xm; i++) { 343 uc = u[j][i].u; 344 vc = u[j][i].v; 345 346 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 347 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 348 349 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 350 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 351 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/ 352 353 stencil[0].i = i; stencil[0].c = 0; entries[0] = -appctx->D1*sy; 354 stencil[1].i = i; stencil[1].c = 0; entries[1] = -appctx->D1*sy; 355 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = -appctx->D1*sx; 356 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = -appctx->D1*sx; 357 stencil[4].i = i; stencil[4].c = 0; entries[4] = 2.0*appctx->D1*(sx + sy) + vc*vc + appctx->gamma + a; 358 stencil[5].i = i; stencil[5].c = 1; entries[5] = 2.0*uc*vc; 359 rowstencil.i = i; rowstencil.c = 0; 360 361 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 362 if (appctx->aijpc) { 363 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 364 } 365 stencil[0].c = 1; entries[0] = -appctx->D2*sy; 366 stencil[1].c = 1; entries[1] = -appctx->D2*sy; 367 stencil[2].c = 1; entries[2] = -appctx->D2*sx; 368 stencil[3].c = 1; entries[3] = -appctx->D2*sx; 369 stencil[4].c = 1; entries[4] = 2.0*appctx->D2*(sx + sy) - 2.0*uc*vc + appctx->gamma + appctx->kappa + a; 370 stencil[5].c = 0; entries[5] = -vc*vc; 371 rowstencil.c = 1; 372 373 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 374 if (appctx->aijpc) { 375 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 376 } 377 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 378 } 379 } 380 381 /* 382 Restore vectors 383 */ 384 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 385 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 386 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 387 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 388 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 389 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 390 if (appctx->aijpc) { 391 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 392 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 393 ierr = MatSetOption(BB,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 394 } 395 PetscFunctionReturn(0); 396 } 397