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