Actual source code: ex15.c
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
8: Boundary conditions:
9: Drichlet BC:
10: At x=0, x=1, y=0, y=1: u = 0.0
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
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
20: */
22: #include <petscdm.h>
23: #include <petscdmda.h>
24: #include <petscts.h>
26: /*
27: User-defined data structures and routines
28: */
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;
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 *);
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; */
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;
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));
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;
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Extract global vectors from DMDA;
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: PetscCall(DMCreateGlobalVector(da, &u));
86: PetscCall(VecDuplicate(u, &r));
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));
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Set initial conditions
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: PetscCall(FormInitialSolution(u, &user));
103: PetscCall(TSSetSolution(ts, u));
104: dt = .01;
105: PetscCall(TSSetTimeStep(ts, dt));
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: }
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Sets various TS parameters from user options
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: PetscCall(TSSetFromOptions(ts));
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Solve nonlinear system
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135: PetscCall(TSSolve(ts, u));
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));
147: PetscCall(PetscFinalize());
148: return 0;
149: }
151: /* --------------------------------------------------------------------- */
152: /*
153: FormIFunction = Udot - RHSFunction
154: */
155: PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
156: {
157: AppCtx *user = (AppCtx *)ctx;
158: DM da = (DM)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;
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));
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");
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));
183: /* Get pointers to vector data */
184: PetscCall(DMDAVecGetArrayRead(da, localU, &uarray));
185: PetscCall(DMDAVecGetArray(da, F, &f));
186: PetscCall(DMDAVecGetArray(da, Udot, &udot));
188: /* Get local grid boundaries */
189: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
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: }
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: }
241: /* --------------------------------------------------------------------- */
242: /*
243: FormIJacobian() - Compute IJacobian = dF/dU + a dF/dUdot
244: This routine is not used with option '-use_coloring'
245: */
246: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *ctx)
247: {
248: PetscInt i, j, Mx, My, xs, ys, xm, ym, nc;
249: AppCtx *user = (AppCtx *)ctx;
250: DM da = (DM)user->da;
251: MatStencil col[5], row;
252: PetscScalar vals[5], hx, hy, sx, sy;
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));
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);
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;
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: }
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: }
335: /* ------------------------------------------------------------------- */
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;
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));
348: hx = 1.0 / (PetscReal)(Mx - 1);
349: hy = 1.0 / (PetscReal)(My - 1);
351: /* Get pointers to vector data */
352: PetscCall(DMDAVecGetArray(da, U, &u));
354: /* Get local grid boundaries */
355: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
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: }
368: /* Restore vectors */
369: PetscCall(DMDAVecRestoreArray(da, U, &u));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: /*TEST
375: test:
376: args: -da_grid_x 20 -da_grid_y 20 -boundary 0 -ts_max_steps 10 -ts_monitor
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
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
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
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
398: TEST*/