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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47: PetscInitialize(&argc, &argv, (char *)0, help);
48: PetscOptionsGetBool(NULL, NULL, "-usemonitor", &usemonitor, NULL);
50: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: Create distributed array (DMDA) to manage parallel grid and vectors
52: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53: 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: DMSetFromOptions(da);
55: DMSetUp(da);
56: DMDASetFieldName(da, 0, "u");
57: DMDASetFieldName(da, 1, "v");
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Extract global vectors from DMDA; then duplicate for remaining
61: vectors that are the same types
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: DMCreateGlobalVector(da, &x);
64: VecDuplicate(x, &r);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Create timestepping solver context
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69: TSCreate(PETSC_COMM_WORLD, &ts);
70: TSSetDM(ts, da);
71: TSSetProblemType(ts, TS_NONLINEAR);
72: TSSetRHSFunction(ts, NULL, FormFunction, da);
74: TSSetMaxTime(ts, 1.0);
75: if (usemonitor) TSMonitorSet(ts, MyTSMonitor, 0, 0);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Customize nonlinear solver
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: TSSetType(ts, TSBEULER);
81: TSGetSNES(ts, &ts_snes);
82: if (usemonitor) {
83: PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf);
84: SNESMonitorSet(ts_snes, (PetscErrorCode(*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscErrorCode(*)(void **))PetscViewerAndFormatDestroy);
85: }
86: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
87: Set initial conditions
88: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: FormInitialSolution(da, x);
90: TSSetTimeStep(ts, .0001);
91: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
92: TSSetSolution(ts, x);
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Set runtime options
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: TSSetFromOptions(ts);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Solve nonlinear system
101: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
102: TSSolve(ts, x);
103: TSGetSolveTime(ts, &ftime);
104: TSGetStepNumber(ts, &steps);
105: 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: VecDestroy(&x);
112: VecDestroy(&r);
113: TSDestroy(&ts);
114: DMDestroy(&da);
116: 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;
140: DMGetLocalVector(da, &localX);
141: 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: DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
157: DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
159: /*
160: Get pointers to vector data
161: */
162: DMDAVecGetArrayDOF(da, localX, &x);
163: DMDAVecGetArrayDOF(da, F, &f);
165: /*
166: Get local grid boundaries
167: */
168: 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: DMDAVecRestoreArrayDOF(da, localX, &x);
193: DMDAVecRestoreArrayDOF(da, F, &f);
194: DMRestoreLocalVector(da, &localX);
195: PetscLogFlops(11.0 * ym * xm);
196: return 0;
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;
207: 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: DMDAVecGetArrayDOF(da, U, &u);
217: /*
218: Get local grid boundaries
219: */
220: 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: DMDAVecRestoreArrayDOF(da, U, &u);
244: return 0;
245: }
247: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
248: {
249: PetscReal norm;
250: MPI_Comm comm;
253: VecNorm(v, NORM_2, &norm);
254: PetscObjectGetComm((PetscObject)ts, &comm);
255: if (step > -1) { /* -1 is used to indicate an interpolated value */
256: PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm);
257: }
258: return 0;
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: {
273: SNESMonitorDefaultShort(snes, its, fnorm, vf);
274: return 0;
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*/