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 */
RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void * ptr)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
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,PetscCtx ctx)92 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat BB, PetscCtx 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 */
IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)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
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,PetscCtx ctx)302 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat BB, PetscCtx 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