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: }
104: /* ------------------------------------------------------------------- */
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: /*
175: RHSJacobian - User-provided routine to compute the Jacobian of
176: the nonlinear right-hand-side function of the ODE.
178: Input Parameters:
179: ts - the TS context
180: t - current time
181: U - global input vector
182: dummy - optional user-defined context, as set by TSetRHSJacobian()
184: Output Parameters:
185: J - Jacobian matrix
186: Jpre - optionally different preconditioning matrix
187: str - flag indicating matrix structure
188: */
189: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx)
190: {
191: DM da;
192: DMDALocalInfo info;
193: PetscInt i, j;
194: PetscReal hx, hy, sx, sy;
196: PetscFunctionBeginUser;
197: PetscCall(TSGetDM(ts, &da));
198: PetscCall(DMDAGetLocalInfo(da, &info));
199: hx = 1.0 / (PetscReal)(info.mx - 1);
200: sx = 1.0 / (hx * hx);
201: hy = 1.0 / (PetscReal)(info.my - 1);
202: sy = 1.0 / (hy * hy);
203: for (j = info.ys; j < info.ys + info.ym; j++) {
204: for (i = info.xs; i < info.xs + info.xm; i++) {
205: PetscInt nc = 0;
206: MatStencil row, col[5];
207: PetscScalar val[5];
208: row.i = i;
209: row.j = j;
210: if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) {
211: col[nc].i = i;
212: col[nc].j = j;
213: val[nc++] = 1.0;
214: } else {
215: col[nc].i = i - 1;
216: col[nc].j = j;
217: val[nc++] = sx;
218: col[nc].i = i + 1;
219: col[nc].j = j;
220: val[nc++] = sx;
221: col[nc].i = i;
222: col[nc].j = j - 1;
223: val[nc++] = sy;
224: col[nc].i = i;
225: col[nc].j = j + 1;
226: val[nc++] = sy;
227: col[nc].i = i;
228: col[nc].j = j;
229: val[nc++] = -2 * sx - 2 * sy;
230: }
231: PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
232: }
233: }
234: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
235: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
236: if (J != Jpre) {
237: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
238: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
239: }
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /* ------------------------------------------------------------------- */
244: PetscErrorCode FormInitialSolution(DM da, Vec U, void *ptr)
245: {
246: AppCtx *user = (AppCtx *)ptr;
247: PetscReal c = user->c;
248: PetscInt i, j, xs, ys, xm, ym, Mx, My;
249: PetscScalar **u;
250: PetscReal hx, hy, x, y, r;
252: PetscFunctionBeginUser;
253: 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));
255: hx = 1.0 / (PetscReal)(Mx - 1);
256: hy = 1.0 / (PetscReal)(My - 1);
258: /* Get pointers to vector data */
259: PetscCall(DMDAVecGetArray(da, U, &u));
261: /* Get local grid boundaries */
262: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
264: /* Compute function over the locally owned part of the grid */
265: for (j = ys; j < ys + ym; j++) {
266: y = j * hy;
267: for (i = xs; i < xs + xm; i++) {
268: x = i * hx;
269: r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
270: if (r < .125) u[j][i] = PetscExpReal(c * r * r * r);
271: else u[j][i] = 0.0;
272: }
273: }
275: /* Restore vectors */
276: PetscCall(DMDAVecRestoreArray(da, U, &u));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*TEST
282: test:
283: args: -ts_max_steps 5 -ts_monitor
285: test:
286: suffix: 2
287: args: -ts_max_steps 5 -ts_monitor
289: test:
290: suffix: 3
291: args: -ts_max_steps 5 -snes_fd_color -ts_monitor
293: TEST*/