1 static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
2 /*
3 Solves the equation
4
5 u_tt - \Delta u = 0
6
7 which we split into two first-order systems
8
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
11
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>
24
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 *);
31
main(int argc,char ** argv)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;
42
43 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
44 Initialize program
45 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
46 PetscFunctionBeginUser;
47 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
48 PetscCall(PetscOptionsGetBool(NULL, NULL, "-usemonitor", &usemonitor, NULL));
49
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"));
58
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));
65
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));
73
74 PetscCall(TSSetMaxTime(ts, 1.0));
75 if (usemonitor) PetscCall(TSMonitorSet(ts, MyTSMonitor, 0, 0));
76
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, (PetscCtxDestroyFn *)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));
93
94 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95 Set runtime options
96 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97 PetscCall(TSSetFromOptions(ts));
98
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"));
106
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));
115
116 PetscCall(PetscFinalize());
117 return 0;
118 }
119 /* ------------------------------------------------------------------- */
120 /*
121 FormFunction - Evaluates nonlinear function, F(x).
122
123 Input Parameters:
124 . ts - the TS context
125 . X - input vector
126 . ptr - optional user-defined context, as set by SNESSetFunction()
127
128 Output Parameter:
129 . F - function vector
130 */
FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void * ptr)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;
138
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));
142
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;*/
149
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));
158
159 /*
160 Get pointers to vector data
161 */
162 PetscCall(DMDAVecGetArrayDOF(da, localX, &x));
163 PetscCall(DMDAVecGetArrayDOF(da, F, &f));
164
165 /*
166 Get local grid boundaries
167 */
168 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
169
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 }
188
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 }
198
199 /* ------------------------------------------------------------------- */
FormInitialSolution(DM da,Vec U)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;
205
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));
208
209 hx = 1.0 / (PetscReal)(Mx - 1);
210 hy = 1.0 / (PetscReal)(My - 1);
211
212 /*
213 Get pointers to vector data
214 */
215 PetscCall(DMDAVecGetArrayDOF(da, U, &u));
216
217 /*
218 Get local grid boundaries
219 */
220 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
221
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 }
239
240 /*
241 Restore vectors
242 */
243 PetscCall(DMDAVecRestoreArrayDOF(da, U, &u));
244 PetscFunctionReturn(PETSC_SUCCESS);
245 }
246
MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscCtx ctx)247 PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscCtx ctx)
248 {
249 PetscReal norm;
250 MPI_Comm comm;
251
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 }
260
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 */
MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat * vf)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
277
278 test:
279 args: -da_grid_x 20 -ts_max_time 3 -ts_time_step 1e-1 -ts_theta_initial_guess_extrapolate 0 -ts_monitor -ksp_monitor_short
280 requires: !single
281
282 test:
283 suffix: 2
284 args: -da_grid_x 20 -ts_max_time 0.11 -ts_time_step 1e-1 -ts_type glle -ts_monitor -ksp_monitor_short
285 requires: !single
286
287 test:
288 suffix: glvis_da_2d_vect
289 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
290 requires: !single
291
292 test:
293 suffix: glvis_da_2d_vect_ll
294 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
295 requires: !single
296
297 TEST*/
298