1 static char help[] = "Time-dependent PDE in 2d. Modified from ex13.c for illustrating how to solve DAEs. \n";
2 /*
3 u_t = uxx + uyy
4 0 < x < 1, 0 < y < 1;
5 At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125
6 u(x,y) = 0.0 if r >= .125
7
8 Boundary conditions:
9 Drichlet BC:
10 At x=0, x=1, y=0, y=1: u = 0.0
11
12 Neumann BC:
13 At x=0, x=1: du(x,y,t)/dx = 0
14 At y=0, y=1: du(x,y,t)/dy = 0
15
16 mpiexec -n 2 ./ex15 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
17 ./ex15 -da_grid_x 40 -da_grid_y 40 -draw_pause .1 -boundary 1 -ts_monitor_draw_solution
18 ./ex15 -da_grid_x 40 -da_grid_y 40 -draw_pause .1 -boundary 1 -Jtype 2 -nstencilpts 9
19
20 */
21
22 #include <petscdm.h>
23 #include <petscdmda.h>
24 #include <petscts.h>
25
26 /*
27 User-defined data structures and routines
28 */
29
30 /* AppCtx: used by FormIFunction() and FormIJacobian() */
31 typedef struct {
32 DM da;
33 PetscInt nstencilpts; /* number of stencil points: 5 or 9 */
34 PetscReal c;
35 PetscInt boundary; /* Type of boundary condition */
36 PetscBool viewJacobian;
37 } AppCtx;
38
39 extern PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
40 extern PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
41 extern PetscErrorCode FormInitialSolution(Vec, void *);
42
main(int argc,char ** argv)43 int main(int argc, char **argv)
44 {
45 TS ts; /* nonlinear solver */
46 Vec u, r; /* solution, residual vectors */
47 Mat J, Jmf = NULL; /* Jacobian matrices */
48 DM da;
49 PetscReal dt;
50 AppCtx user; /* user-defined work context */
51 SNES snes;
52 PetscInt Jtype; /* Jacobian type
53 0: user provide Jacobian;
54 1: slow finite difference;
55 2: fd with coloring; */
56
57 PetscFunctionBeginUser;
58 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
59 /* Initialize user application context */
60 user.da = NULL;
61 user.nstencilpts = 5;
62 user.c = -30.0;
63 user.boundary = 0; /* 0: Drichlet BC; 1: Neumann BC */
64 user.viewJacobian = PETSC_FALSE;
65
66 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nstencilpts", &user.nstencilpts, NULL));
67 PetscCall(PetscOptionsGetInt(NULL, NULL, "-boundary", &user.boundary, NULL));
68 PetscCall(PetscOptionsHasName(NULL, NULL, "-viewJacobian", &user.viewJacobian));
69
70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 Create distributed array (DMDA) to manage parallel grid and vectors
72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73 if (user.nstencilpts == 5) {
74 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 11, 11, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
75 } else if (user.nstencilpts == 9) {
76 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 11, 11, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
77 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "nstencilpts %" PetscInt_FMT " is not supported", user.nstencilpts);
78 PetscCall(DMSetFromOptions(da));
79 PetscCall(DMSetUp(da));
80 user.da = da;
81
82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 Extract global vectors from DMDA;
84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85 PetscCall(DMCreateGlobalVector(da, &u));
86 PetscCall(VecDuplicate(u, &r));
87
88 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89 Create timestepping solver context
90 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
92 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
93 PetscCall(TSSetType(ts, TSBEULER));
94 PetscCall(TSSetDM(ts, da));
95 PetscCall(TSSetIFunction(ts, r, FormIFunction, &user));
96 PetscCall(TSSetMaxTime(ts, 1.0));
97 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
98
99 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100 Set initial conditions
101 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102 PetscCall(FormInitialSolution(u, &user));
103 PetscCall(TSSetSolution(ts, u));
104 dt = .01;
105 PetscCall(TSSetTimeStep(ts, dt));
106
107 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108 Set Jacobian evaluation routine
109 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110 PetscCall(DMSetMatType(da, MATAIJ));
111 PetscCall(DMCreateMatrix(da, &J));
112 Jtype = 0;
113 PetscCall(PetscOptionsGetInt(NULL, NULL, "-Jtype", &Jtype, NULL));
114 if (Jtype == 0) { /* use user provided Jacobian evaluation routine */
115 PetscCheck(user.nstencilpts == 5, PETSC_COMM_WORLD, PETSC_ERR_SUP, "user Jacobian routine FormIJacobian() does not support nstencilpts=%" PetscInt_FMT, user.nstencilpts);
116 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));
117 } else { /* use finite difference Jacobian J as preconditioner and '-snes_mf_operator' for Mat*vec */
118 PetscCall(TSGetSNES(ts, &snes));
119 PetscCall(MatCreateSNESMF(snes, &Jmf));
120 if (Jtype == 1) { /* slow finite difference J; */
121 PetscCall(SNESSetJacobian(snes, Jmf, J, SNESComputeJacobianDefault, NULL));
122 } else if (Jtype == 2) { /* Use coloring to compute finite difference J efficiently */
123 PetscCall(SNESSetJacobian(snes, Jmf, J, SNESComputeJacobianDefaultColor, 0));
124 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Jtype is not supported");
125 }
126
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 Sets various TS parameters from user options
129 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130 PetscCall(TSSetFromOptions(ts));
131
132 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133 Solve nonlinear system
134 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135 PetscCall(TSSolve(ts, u));
136
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 Free work space.
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140 PetscCall(MatDestroy(&J));
141 PetscCall(MatDestroy(&Jmf));
142 PetscCall(VecDestroy(&u));
143 PetscCall(VecDestroy(&r));
144 PetscCall(TSDestroy(&ts));
145 PetscCall(DMDestroy(&da));
146
147 PetscCall(PetscFinalize());
148 return 0;
149 }
150
151 /* --------------------------------------------------------------------- */
152 /*
153 FormIFunction = Udot - RHSFunction
154 */
FormIFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)155 PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
156 {
157 AppCtx *user = (AppCtx *)ctx;
158 DM da = user->da;
159 PetscInt i, j, Mx, My, xs, ys, xm, ym;
160 PetscReal hx, hy, sx, sy;
161 PetscScalar u, uxx, uyy, **uarray, **f, **udot;
162 Vec localU;
163
164 PetscFunctionBeginUser;
165 PetscCall(DMGetLocalVector(da, &localU));
166 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));
167
168 hx = 1.0 / (PetscReal)(Mx - 1);
169 sx = 1.0 / (hx * hx);
170 hy = 1.0 / (PetscReal)(My - 1);
171 sy = 1.0 / (hy * hy);
172 PetscCheck(user->nstencilpts != 9 || hx == hy, PETSC_COMM_WORLD, PETSC_ERR_SUP, "hx must equal hy when nstencilpts = 9 for this example");
173
174 /*
175 Scatter ghost points to local vector,using the 2-step process
176 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
177 By placing code between these two statements, computations can be
178 done while messages are in transition.
179 */
180 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
181 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
182
183 /* Get pointers to vector data */
184 PetscCall(DMDAVecGetArrayRead(da, localU, &uarray));
185 PetscCall(DMDAVecGetArray(da, F, &f));
186 PetscCall(DMDAVecGetArray(da, Udot, &udot));
187
188 /* Get local grid boundaries */
189 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
190
191 /* Compute function over the locally owned part of the grid */
192 for (j = ys; j < ys + ym; j++) {
193 for (i = xs; i < xs + xm; i++) {
194 /* Boundary conditions */
195 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
196 if (user->boundary == 0) { /* Drichlet BC */
197 f[j][i] = uarray[j][i]; /* F = U */
198 } else { /* Neumann BC */
199 if (i == 0 && j == 0) { /* SW corner */
200 f[j][i] = uarray[j][i] - uarray[j + 1][i + 1];
201 } else if (i == Mx - 1 && j == 0) { /* SE corner */
202 f[j][i] = uarray[j][i] - uarray[j + 1][i - 1];
203 } else if (i == 0 && j == My - 1) { /* NW corner */
204 f[j][i] = uarray[j][i] - uarray[j - 1][i + 1];
205 } else if (i == Mx - 1 && j == My - 1) { /* NE corner */
206 f[j][i] = uarray[j][i] - uarray[j - 1][i - 1];
207 } else if (i == 0) { /* Left */
208 f[j][i] = uarray[j][i] - uarray[j][i + 1];
209 } else if (i == Mx - 1) { /* Right */
210 f[j][i] = uarray[j][i] - uarray[j][i - 1];
211 } else if (j == 0) { /* Bottom */
212 f[j][i] = uarray[j][i] - uarray[j + 1][i];
213 } else if (j == My - 1) { /* Top */
214 f[j][i] = uarray[j][i] - uarray[j - 1][i];
215 }
216 }
217 } else { /* Interior */
218 u = uarray[j][i];
219 /* 5-point stencil */
220 uxx = (-2.0 * u + uarray[j][i - 1] + uarray[j][i + 1]);
221 uyy = (-2.0 * u + uarray[j - 1][i] + uarray[j + 1][i]);
222 if (user->nstencilpts == 9) {
223 /* 9-point stencil: assume hx=hy */
224 uxx = 2.0 * uxx / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] + uarray[j + 1][i - 1] + uarray[j + 1][i + 1]) - 2.0 * u) / 6.0;
225 uyy = 2.0 * uyy / 3.0 + (0.5 * (uarray[j - 1][i - 1] + uarray[j - 1][i + 1] + uarray[j + 1][i - 1] + uarray[j + 1][i + 1]) - 2.0 * u) / 6.0;
226 }
227 f[j][i] = udot[j][i] - (uxx * sx + uyy * sy);
228 }
229 }
230 }
231
232 /* Restore vectors */
233 PetscCall(DMDAVecRestoreArrayRead(da, localU, &uarray));
234 PetscCall(DMDAVecRestoreArray(da, F, &f));
235 PetscCall(DMDAVecRestoreArray(da, Udot, &udot));
236 PetscCall(DMRestoreLocalVector(da, &localU));
237 PetscCall(PetscLogFlops(11.0 * ym * xm));
238 PetscFunctionReturn(PETSC_SUCCESS);
239 }
240
241 /* --------------------------------------------------------------------- */
242 /*
243 FormIJacobian() - Compute IJacobian = dF/dU + a dF/dUdot
244 This routine is not used with option '-use_coloring'
245 */
FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,PetscCtx ctx)246 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, PetscCtx ctx)
247 {
248 PetscInt i, j, Mx, My, xs, ys, xm, ym, nc;
249 AppCtx *user = (AppCtx *)ctx;
250 DM da = user->da;
251 MatStencil col[5], row;
252 PetscScalar vals[5], hx, hy, sx, sy;
253
254 PetscFunctionBeginUser;
255 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));
256 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
257
258 hx = 1.0 / (PetscReal)(Mx - 1);
259 sx = 1.0 / (hx * hx);
260 hy = 1.0 / (PetscReal)(My - 1);
261 sy = 1.0 / (hy * hy);
262
263 for (j = ys; j < ys + ym; j++) {
264 for (i = xs; i < xs + xm; i++) {
265 nc = 0;
266 row.j = j;
267 row.i = i;
268 if (user->boundary == 0 && (i == 0 || i == Mx - 1 || j == 0 || j == My - 1)) {
269 col[nc].j = j;
270 col[nc].i = i;
271 vals[nc++] = 1.0;
272
273 } else if (user->boundary > 0 && i == 0) { /* Left Neumann */
274 col[nc].j = j;
275 col[nc].i = i;
276 vals[nc++] = 1.0;
277 col[nc].j = j;
278 col[nc].i = i + 1;
279 vals[nc++] = -1.0;
280 } else if (user->boundary > 0 && i == Mx - 1) { /* Right Neumann */
281 col[nc].j = j;
282 col[nc].i = i;
283 vals[nc++] = 1.0;
284 col[nc].j = j;
285 col[nc].i = i - 1;
286 vals[nc++] = -1.0;
287 } else if (user->boundary > 0 && j == 0) { /* Bottom Neumann */
288 col[nc].j = j;
289 col[nc].i = i;
290 vals[nc++] = 1.0;
291 col[nc].j = j + 1;
292 col[nc].i = i;
293 vals[nc++] = -1.0;
294 } else if (user->boundary > 0 && j == My - 1) { /* Top Neumann */
295 col[nc].j = j;
296 col[nc].i = i;
297 vals[nc++] = 1.0;
298 col[nc].j = j - 1;
299 col[nc].i = i;
300 vals[nc++] = -1.0;
301 } else { /* Interior */
302 col[nc].j = j - 1;
303 col[nc].i = i;
304 vals[nc++] = -sy;
305 col[nc].j = j;
306 col[nc].i = i - 1;
307 vals[nc++] = -sx;
308 col[nc].j = j;
309 col[nc].i = i;
310 vals[nc++] = 2.0 * (sx + sy) + a;
311 col[nc].j = j;
312 col[nc].i = i + 1;
313 vals[nc++] = -sx;
314 col[nc].j = j + 1;
315 col[nc].i = i;
316 vals[nc++] = -sy;
317 }
318 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, vals, INSERT_VALUES));
319 }
320 }
321 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
322 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
323 if (J != Jpre) {
324 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
325 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
326 }
327
328 if (user->viewJacobian) {
329 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)Jpre), "Jpre:\n"));
330 PetscCall(MatView(Jpre, PETSC_VIEWER_STDOUT_WORLD));
331 }
332 PetscFunctionReturn(PETSC_SUCCESS);
333 }
334
335 /* ------------------------------------------------------------------- */
FormInitialSolution(Vec U,void * ptr)336 PetscErrorCode FormInitialSolution(Vec U, void *ptr)
337 {
338 AppCtx *user = (AppCtx *)ptr;
339 DM da = user->da;
340 PetscReal c = user->c;
341 PetscInt i, j, xs, ys, xm, ym, Mx, My;
342 PetscScalar **u;
343 PetscReal hx, hy, x, y, r;
344
345 PetscFunctionBeginUser;
346 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));
347
348 hx = 1.0 / (PetscReal)(Mx - 1);
349 hy = 1.0 / (PetscReal)(My - 1);
350
351 /* Get pointers to vector data */
352 PetscCall(DMDAVecGetArray(da, U, &u));
353
354 /* Get local grid boundaries */
355 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
356
357 /* Compute function over the locally owned part of the grid */
358 for (j = ys; j < ys + ym; j++) {
359 y = j * hy;
360 for (i = xs; i < xs + xm; i++) {
361 x = i * hx;
362 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
363 if (r < .125) u[j][i] = PetscExpReal(c * r * r * r);
364 else u[j][i] = 0.0;
365 }
366 }
367
368 /* Restore vectors */
369 PetscCall(DMDAVecRestoreArray(da, U, &u));
370 PetscFunctionReturn(PETSC_SUCCESS);
371 }
372
373 /*TEST
374
375 test:
376 args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -ts_monitor
377
378 test:
379 suffix: 2
380 args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -Jtype 2 -ts_monitor
381
382 test:
383 suffix: 3
384 requires: !single
385 args: -da_grid_x 20 -da_grid_y 20 -boundary 1 -ts_max_steps 10 -ts_monitor
386
387 test:
388 suffix: 4
389 requires: !single
390 nsize: 2
391 args: -da_grid_x 20 -da_grid_y 20 -boundary 1 -ts_max_steps 10 -ts_monitor
392
393 test:
394 suffix: 5
395 nsize: 1
396 args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -Jtype 1 -ts_monitor
397
398 TEST*/
399