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