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