Actual source code: ex13.c
1: static char help[] = "Time-dependent PDE in 2d. Simplified from ex7.c for illustrating how to use TS on a structured domain. \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: mpiexec -n 2 ./ex13 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
9: mpiexec -n 1 ./ex13 -snes_fd_color -ts_monitor_draw_solution
10: mpiexec -n 2 ./ex13 -ts_type sundials -ts_monitor
11: */
13: #include <petscdm.h>
14: #include <petscdmda.h>
15: #include <petscts.h>
17: /*
18: User-defined data structures and routines
19: */
20: typedef struct {
21: PetscReal c;
22: } AppCtx;
24: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
25: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
26: extern PetscErrorCode FormInitialSolution(DM, Vec, void *);
28: int main(int argc, char **argv)
29: {
30: TS ts; /* nonlinear solver */
31: Vec u, r; /* solution, residual vector */
32: Mat J; /* Jacobian matrix */
33: PetscInt steps; /* iterations for convergence */
34: DM da;
35: PetscReal ftime, dt;
36: AppCtx user; /* user-defined work context */
38: PetscFunctionBeginUser;
39: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
40: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
41: Create distributed array (DMDA) to manage parallel grid and vectors
42: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 8, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
44: PetscCall(DMSetFromOptions(da));
45: PetscCall(DMSetUp(da));
47: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48: Extract global vectors from DMDA;
49: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50: PetscCall(DMCreateGlobalVector(da, &u));
51: PetscCall(VecDuplicate(u, &r));
53: /* Initialize user application context */
54: user.c = -30.0;
56: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: Create timestepping solver context
58: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
60: PetscCall(TSSetDM(ts, da));
61: PetscCall(TSSetType(ts, TSBEULER));
62: PetscCall(TSSetRHSFunction(ts, r, RHSFunction, &user));
64: /* Set Jacobian */
65: PetscCall(DMSetMatType(da, MATAIJ));
66: PetscCall(DMCreateMatrix(da, &J));
67: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, NULL));
69: ftime = 1.0;
70: PetscCall(TSSetMaxTime(ts, ftime));
71: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Set initial conditions
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: PetscCall(FormInitialSolution(da, u, &user));
77: dt = .01;
78: PetscCall(TSSetTimeStep(ts, dt));
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Set runtime options
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83: PetscCall(TSSetFromOptions(ts));
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Solve nonlinear system
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: PetscCall(TSSolve(ts, u));
89: PetscCall(TSGetSolveTime(ts, &ftime));
90: PetscCall(TSGetStepNumber(ts, &steps));
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Free work space.
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95: PetscCall(MatDestroy(&J));
96: PetscCall(VecDestroy(&u));
97: PetscCall(VecDestroy(&r));
98: PetscCall(TSDestroy(&ts));
99: PetscCall(DMDestroy(&da));
101: PetscCall(PetscFinalize());
102: return 0;
103: }
105: /*
106: RHSFunction - Evaluates nonlinear function, F(u).
108: Input Parameters:
109: . ts - the TS context
110: . U - input vector
111: . ptr - optional user-defined context, as set by TSSetFunction()
113: Output Parameter:
114: . F - function vector
115: */
116: PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec U, Vec F, void *ptr)
117: {
118: /* PETSC_UNUSED AppCtx *user=(AppCtx*)ptr; */
119: DM da;
120: PetscInt i, j, Mx, My, xs, ys, xm, ym;
121: PetscReal two = 2.0, hx, hy, sx, sy;
122: PetscScalar u, uxx, uyy, **uarray, **f;
123: Vec localU;
125: PetscFunctionBeginUser;
126: PetscCall(TSGetDM(ts, &da));
127: PetscCall(DMGetLocalVector(da, &localU));
128: 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));
130: hx = 1.0 / (PetscReal)(Mx - 1);
131: sx = 1.0 / (hx * hx);
132: hy = 1.0 / (PetscReal)(My - 1);
133: sy = 1.0 / (hy * hy);
135: /*
136: Scatter ghost points to local vector,using the 2-step process
137: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
138: By placing code between these two statements, computations can be
139: done while messages are in transition.
140: */
141: PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
142: PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
144: /* Get pointers to vector data */
145: PetscCall(DMDAVecGetArrayRead(da, localU, &uarray));
146: PetscCall(DMDAVecGetArray(da, F, &f));
148: /* Get local grid boundaries */
149: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
151: /* Compute function over the locally owned part of the grid */
152: for (j = ys; j < ys + ym; j++) {
153: for (i = xs; i < xs + xm; i++) {
154: if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
155: f[j][i] = uarray[j][i];
156: continue;
157: }
158: u = uarray[j][i];
159: uxx = (-two * u + uarray[j][i - 1] + uarray[j][i + 1]) * sx;
160: uyy = (-two * u + uarray[j - 1][i] + uarray[j + 1][i]) * sy;
161: f[j][i] = uxx + uyy;
162: }
163: }
165: /* Restore vectors */
166: PetscCall(DMDAVecRestoreArrayRead(da, localU, &uarray));
167: PetscCall(DMDAVecRestoreArray(da, F, &f));
168: PetscCall(DMRestoreLocalVector(da, &localU));
169: PetscCall(PetscLogFlops(11.0 * ym * xm));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*
174: RHSJacobian - User-provided routine to compute the Jacobian of
175: the nonlinear right-hand-side function of the ODE.
177: Input Parameters:
178: ts - the TS context
179: t - current time
180: U - global input vector
181: dummy - optional user-defined context, as set by TSetRHSJacobian()
183: Output Parameters:
184: J - Jacobian matrix
185: Jpre - optionally different preconditioning matrix
187: */
188: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx)
189: {
190: DM da;
191: DMDALocalInfo info;
192: PetscInt i, j;
193: PetscReal hx, hy, sx, sy;
195: PetscFunctionBeginUser;
196: PetscCall(TSGetDM(ts, &da));
197: PetscCall(DMDAGetLocalInfo(da, &info));
198: hx = 1.0 / (PetscReal)(info.mx - 1);
199: sx = 1.0 / (hx * hx);
200: hy = 1.0 / (PetscReal)(info.my - 1);
201: sy = 1.0 / (hy * hy);
202: for (j = info.ys; j < info.ys + info.ym; j++) {
203: for (i = info.xs; i < info.xs + info.xm; i++) {
204: PetscInt nc = 0;
205: MatStencil row, col[5];
206: PetscScalar val[5];
207: row.i = i;
208: row.j = j;
209: if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) {
210: col[nc].i = i;
211: col[nc].j = j;
212: val[nc++] = 1.0;
213: } else {
214: col[nc].i = i - 1;
215: col[nc].j = j;
216: val[nc++] = sx;
217: col[nc].i = i + 1;
218: col[nc].j = j;
219: val[nc++] = sx;
220: col[nc].i = i;
221: col[nc].j = j - 1;
222: val[nc++] = sy;
223: col[nc].i = i;
224: col[nc].j = j + 1;
225: val[nc++] = sy;
226: col[nc].i = i;
227: col[nc].j = j;
228: val[nc++] = -2 * sx - 2 * sy;
229: }
230: PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
231: }
232: }
233: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
234: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
235: if (J != Jpre) {
236: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
237: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
238: }
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: PetscErrorCode FormInitialSolution(DM da, Vec U, void *ptr)
243: {
244: AppCtx *user = (AppCtx *)ptr;
245: PetscReal c = user->c;
246: PetscInt i, j, xs, ys, xm, ym, Mx, My;
247: PetscScalar **u;
248: PetscReal hx, hy, x, y, r;
250: PetscFunctionBeginUser;
251: 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));
253: hx = 1.0 / (PetscReal)(Mx - 1);
254: hy = 1.0 / (PetscReal)(My - 1);
256: /* Get pointers to vector data */
257: PetscCall(DMDAVecGetArray(da, U, &u));
259: /* Get local grid boundaries */
260: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
262: /* Compute function over the locally owned part of the grid */
263: for (j = ys; j < ys + ym; j++) {
264: y = j * hy;
265: for (i = xs; i < xs + xm; i++) {
266: x = i * hx;
267: r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
268: if (r < .125) u[j][i] = PetscExpReal(c * r * r * r);
269: else u[j][i] = 0.0;
270: }
271: }
273: /* Restore vectors */
274: PetscCall(DMDAVecRestoreArray(da, U, &u));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: /*TEST
280: test:
281: args: -ts_max_steps 5 -ts_monitor
283: test:
284: suffix: 2
285: args: -ts_max_steps 5 -ts_monitor
287: test:
288: suffix: 3
289: args: -ts_max_steps 5 -snes_fd_color -ts_monitor
291: TEST*/