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