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