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