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 */
40: PetscInitialize(&argc, &argv, (char *)0, help);
41: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: Create distributed array (DMDA) to manage parallel grid and vectors
43: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
44: 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: DMSetFromOptions(da);
46: DMSetUp(da);
48: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49: Extract global vectors from DMDA;
50: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51: DMCreateGlobalVector(da, &u);
52: VecDuplicate(u, &r);
54: /* Initialize user application context */
55: user.c = -30.0;
57: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: Create timestepping solver context
59: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60: TSCreate(PETSC_COMM_WORLD, &ts);
61: TSSetDM(ts, da);
62: TSSetType(ts, TSBEULER);
63: TSSetRHSFunction(ts, r, RHSFunction, &user);
65: /* Set Jacobian */
66: DMSetMatType(da, MATAIJ);
67: DMCreateMatrix(da, &J);
68: TSSetRHSJacobian(ts, J, J, RHSJacobian, NULL);
70: ftime = 1.0;
71: TSSetMaxTime(ts, ftime);
72: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Set initial conditions
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: FormInitialSolution(da, u, &user);
78: dt = .01;
79: TSSetTimeStep(ts, dt);
81: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82: Set runtime options
83: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: TSSetFromOptions(ts);
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Solve nonlinear system
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: TSSolve(ts, u);
90: TSGetSolveTime(ts, &ftime);
91: TSGetStepNumber(ts, &steps);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Free work space.
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: MatDestroy(&J);
97: VecDestroy(&u);
98: VecDestroy(&r);
99: TSDestroy(&ts);
100: DMDestroy(&da);
102: 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;
127: TSGetDM(ts, &da);
128: DMGetLocalVector(da, &localU);
129: 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: DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU);
143: DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU);
145: /* Get pointers to vector data */
146: DMDAVecGetArrayRead(da, localU, &uarray);
147: DMDAVecGetArray(da, F, &f);
149: /* Get local grid boundaries */
150: 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: DMDAVecRestoreArrayRead(da, localU, &uarray);
168: DMDAVecRestoreArray(da, F, &f);
169: DMRestoreLocalVector(da, &localU);
170: PetscLogFlops(11.0 * ym * xm);
171: return 0;
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;
198: TSGetDM(ts, &da);
199: 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: MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES);
233: }
234: }
235: MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
236: MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
237: if (J != Jpre) {
238: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
239: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
240: }
241: return 0;
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;
254: 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: DMDAVecGetArray(da, U, &u);
262: /* Get local grid boundaries */
263: 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: DMDAVecRestoreArray(da, U, &u);
278: return 0;
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*/