Actual source code: ex12.c
1: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
2: /*
3: Solves the equation
5: u_tt - \Delta u = 0
7: which we split into two first-order systems
9: u_t - v = 0 so that F(u,v).u = v
10: v_t - \Delta u = 0 so that F(u,v).v = Delta u
12: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
13: Include "petscts.h" so that we can use the PETSc ODE integrators. Note that this
14: file automatically includes:
15: petscsys.h - base PETSc routines
16: petscis.h - index sets
17: petscvec.h - vectors
18: petscmat.h - matrices
19: petscpc.h - preconditioners
20: petscksp.h - linear solvers
21: */
22: #include <petscdm.h>
23: #include <petscdmda.h>
24: #include <petscts.h>
26: /*
27: User-defined routines
28: */
29: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec);
30: extern PetscErrorCode MyTSMonitor(TS, PetscInt, PetscReal, Vec, void *);
31: extern PetscErrorCode MySNESMonitor(SNES, PetscInt, PetscReal, PetscViewerAndFormat *);
33: int main(int argc, char **argv)
34: {
35: TS ts; /* nonlinear solver */
36: Vec x, r; /* solution, residual vectors */
37: PetscInt steps; /* iterations for convergence */
38: DM da;
39: PetscReal ftime;
40: SNES ts_snes;
41: PetscBool usemonitor = PETSC_TRUE;
42: PetscViewerAndFormat *vf;
44: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: Initialize program
46: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: PetscFunctionBeginUser;
48: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
49: PetscCall(PetscOptionsGetBool(NULL, NULL, "-usemonitor", &usemonitor, NULL));
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create distributed array (DMDA) to manage parallel grid and vectors
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 8, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, NULL, &da));
55: PetscCall(DMSetFromOptions(da));
56: PetscCall(DMSetUp(da));
57: PetscCall(DMDASetFieldName(da, 0, "u"));
58: PetscCall(DMDASetFieldName(da, 1, "v"));
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Extract global vectors from DMDA; then duplicate for remaining
62: vectors that are the same types
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: PetscCall(DMCreateGlobalVector(da, &x));
65: PetscCall(VecDuplicate(x, &r));
67: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68: Create timestepping solver context
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
71: PetscCall(TSSetDM(ts, da));
72: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
73: PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, da));
75: PetscCall(TSSetMaxTime(ts, 1.0));
76: if (usemonitor) PetscCall(TSMonitorSet(ts, MyTSMonitor, NULL, NULL));
78: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79: Customize nonlinear solver
80: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: PetscCall(TSSetType(ts, TSBEULER));
82: PetscCall(TSGetSNES(ts, &ts_snes));
83: if (usemonitor) {
84: PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
85: PetscCall(SNESMonitorSet(ts_snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));
86: }
87: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: Set initial conditions
89: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
90: PetscCall(FormInitialSolution(da, x));
91: PetscCall(TSSetTimeStep(ts, .0001));
92: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
93: PetscCall(TSSetSolution(ts, x));
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Set runtime options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscCall(TSSetFromOptions(ts));
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Solve nonlinear system
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: PetscCall(TSSolve(ts, x));
104: PetscCall(TSGetSolveTime(ts, &ftime));
105: PetscCall(TSGetStepNumber(ts, &steps));
106: PetscCall(VecViewFromOptions(x, NULL, "-final_sol"));
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Free work space. All PETSc objects should be destroyed when they
110: are no longer needed.
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: PetscCall(VecDestroy(&x));
113: PetscCall(VecDestroy(&r));
114: PetscCall(TSDestroy(&ts));
115: PetscCall(DMDestroy(&da));
117: PetscCall(PetscFinalize());
118: return 0;
119: }
120: /* ------------------------------------------------------------------- */
121: /*
122: FormFunction - Evaluates nonlinear function, F(x).
124: Input Parameters:
125: . ts - the TS context
126: . X - input vector
127: . ptr - optional user-defined context, as set by SNESSetFunction()
129: Output Parameter:
130: . F - function vector
131: */
132: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
133: {
134: DM da = (DM)ptr;
135: PetscInt i, j, Mx, My, xs, ys, xm, ym;
136: PetscReal hx, hy, /*hxdhy,hydhx,*/ sx, sy;
137: PetscScalar u, uxx, uyy, v, ***x, ***f;
138: Vec localX;
140: PetscFunctionBeginUser;
141: PetscCall(DMGetLocalVector(da, &localX));
142: 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));
144: hx = 1.0 / (PetscReal)(Mx - 1);
145: sx = 1.0 / (hx * hx);
146: hy = 1.0 / (PetscReal)(My - 1);
147: sy = 1.0 / (hy * hy);
148: /*hxdhy = hx/hy;*/
149: /*hydhx = hy/hx;*/
151: /*
152: Scatter ghost points to local vector,using the 2-step process
153: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
154: By placing code between these two statements, computations can be
155: done while messages are in transition.
156: */
157: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
158: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
160: /*
161: Get pointers to vector data
162: */
163: PetscCall(DMDAVecGetArrayDOF(da, localX, &x));
164: PetscCall(DMDAVecGetArrayDOF(da, F, &f));
166: /*
167: Get local grid boundaries
168: */
169: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
171: /*
172: Compute function over the locally owned part of the grid
173: */
174: for (j = ys; j < ys + ym; j++) {
175: for (i = xs; i < xs + xm; i++) {
176: if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
177: f[j][i][0] = x[j][i][0];
178: f[j][i][1] = x[j][i][1];
179: continue;
180: }
181: u = x[j][i][0];
182: v = x[j][i][1];
183: uxx = (-2.0 * u + x[j][i - 1][0] + x[j][i + 1][0]) * sx;
184: uyy = (-2.0 * u + x[j - 1][i][0] + x[j + 1][i][0]) * sy;
185: f[j][i][0] = v;
186: f[j][i][1] = uxx + uyy;
187: }
188: }
190: /*
191: Restore vectors
192: */
193: PetscCall(DMDAVecRestoreArrayDOF(da, localX, &x));
194: PetscCall(DMDAVecRestoreArrayDOF(da, F, &f));
195: PetscCall(DMRestoreLocalVector(da, &localX));
196: PetscCall(PetscLogFlops(11.0 * ym * xm));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: /* ------------------------------------------------------------------- */
201: PetscErrorCode FormInitialSolution(DM da, Vec U)
202: {
203: PetscInt i, j, xs, ys, xm, ym, Mx, My;
204: PetscScalar ***u;
205: PetscReal hx, hy, x, y, r;
207: PetscFunctionBeginUser;
208: 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));
210: hx = 1.0 / (PetscReal)(Mx - 1);
211: hy = 1.0 / (PetscReal)(My - 1);
213: /*
214: Get pointers to vector data
215: */
216: PetscCall(DMDAVecGetArrayDOF(da, U, &u));
218: /*
219: Get local grid boundaries
220: */
221: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
223: /*
224: Compute function over the locally owned part of the grid
225: */
226: for (j = ys; j < ys + ym; j++) {
227: y = j * hy;
228: for (i = xs; i < xs + xm; i++) {
229: x = i * hx;
230: r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
231: if (r < .125) {
232: u[j][i][0] = PetscExpReal(-30.0 * r * r * r);
233: u[j][i][1] = 0.0;
234: } else {
235: u[j][i][0] = 0.0;
236: u[j][i][1] = 0.0;
237: }
238: }
239: }
241: /*
242: Restore vectors
243: */
244: PetscCall(DMDAVecRestoreArrayDOF(da, U, &u));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscCtx ctx)
249: {
250: PetscReal norm;
251: MPI_Comm comm;
253: PetscFunctionBeginUser;
254: PetscCall(VecNorm(v, NORM_2, &norm));
255: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
256: if (step > -1) { /* -1 is used to indicate an interpolated value */
257: PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
258: }
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*
263: MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
264: Input Parameters:
265: snes - the SNES context
266: its - iteration number
267: fnorm - 2-norm function value (may be estimated)
268: ctx - optional user-defined context for private data for the
269: monitor routine, as set by SNESMonitorSet()
270: */
271: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
272: {
273: PetscFunctionBeginUser;
274: PetscCall(SNESMonitorDefault(snes, its, fnorm, vf));
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
277: /*TEST
279: test:
280: args: -da_grid_x 20 -ts_max_time 3 -ts_time_step 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor
281: requires: !single
283: test:
284: suffix: 2
285: args: -da_grid_x 20 -ts_max_time 0.11 -ts_time_step 1e-1 -ts_type glle -ts_monitor -ksp_monitor
286: requires: !single
288: test:
289: suffix: glvis_da_2d_vect
290: args: -usemonitor 0 -da_grid_x 20 -ts_max_time 0.3 -ts_time_step 1e-1 -ts_type glle -final_sol glvis: -viewer_glvis_dm_da_bs 2,0
291: requires: !single
293: test:
294: suffix: glvis_da_2d_vect_ll
295: args: -usemonitor 0 -da_grid_x 20 -ts_max_time 0.3 -ts_time_step 1e-1 -ts_type glle -final_sol glvis: -viewer_glvis_dm_da_bs 2,0 -viewer_glvis_dm_da_ll
296: requires: !single
298: TEST*/