1 static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";
2
3 /*
4 Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
5 Include "petscts.h" so that we can use SNES solvers. Note that this
6 file automatically includes:
7 petscsys.h - base PETSc routines petscvec.h - vectors
8 petscmat.h - matrices
9 petscis.h - index sets petscksp.h - Krylov subspace methods
10 petscviewer.h - viewers petscpc.h - preconditioners
11 petscksp.h - linear solvers
12 */
13 #include <petscdm.h>
14 #include <petscdmda.h>
15 #include <petscts.h>
16
17 /*
18 User-defined routines
19 */
20 extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec);
21 extern PetscErrorCode MyTSMonitor(TS, PetscInt, PetscReal, Vec, void *);
22 extern PetscErrorCode MySNESMonitor(SNES, PetscInt, PetscReal, PetscViewerAndFormat *);
23
main(int argc,char ** argv)24 int main(int argc, char **argv)
25 {
26 TS ts; /* time integrator */
27 SNES snes;
28 Vec x, r; /* solution, residual vectors */
29 DM da;
30 PetscViewerAndFormat *vf;
31
32 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
33 Initialize program
34 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
35 PetscFunctionBeginUser;
36 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
37 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38 Create distributed array (DMDA) to manage parallel grid and vectors
39 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 8, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
41 PetscCall(DMSetFromOptions(da));
42 PetscCall(DMSetUp(da));
43
44 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45 Extract global vectors from DMDA; then duplicate for remaining
46 vectors that are the same types
47 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
48 PetscCall(DMCreateGlobalVector(da, &x));
49 PetscCall(VecDuplicate(x, &r));
50
51 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52 Create timestepping solver context
53 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
55 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
56 PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, da));
57
58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59 Create matrix data structure; set Jacobian evaluation routine
60
61 Set Jacobian matrix data structure and default Jacobian evaluation
62 routine. User can override with:
63 -snes_mf : matrix-free Newton-Krylov method with no preconditioning
64 (unless user explicitly sets preconditioner)
65 -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
66 but use matrix-free approx for Jacobian-vector
67 products within Newton-Krylov method
68
69 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70
71 PetscCall(TSSetMaxTime(ts, 1.0));
72 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
73 PetscCall(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL));
74 PetscCall(TSSetDM(ts, da));
75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76 Customize nonlinear solver
77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78 PetscCall(TSSetType(ts, TSBEULER));
79 PetscCall(TSGetSNES(ts, &snes));
80 PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
81 PetscCall(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));
82
83 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84 Set initial conditions
85 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86 PetscCall(FormInitialSolution(da, x));
87 PetscCall(TSSetTimeStep(ts, .0001));
88 PetscCall(TSSetSolution(ts, x));
89
90 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91 Set runtime options
92 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93 PetscCall(TSSetFromOptions(ts));
94
95 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96 Solve nonlinear system
97 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98 PetscCall(TSSolve(ts, x));
99
100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101 Free work space. All PETSc objects should be destroyed when they
102 are no longer needed.
103 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104 PetscCall(VecDestroy(&x));
105 PetscCall(VecDestroy(&r));
106 PetscCall(TSDestroy(&ts));
107 PetscCall(DMDestroy(&da));
108
109 PetscCall(PetscFinalize());
110 return 0;
111 }
112 /* ------------------------------------------------------------------- */
113 /*
114 FormFunction - Evaluates nonlinear function, F(x).
115
116 Input Parameters:
117 . ts - the TS context
118 . X - input vector
119 . ptr - optional user-defined context, as set by SNESSetFunction()
120
121 Output Parameter:
122 . F - function vector
123 */
FormFunction(TS ts,PetscReal ftime,Vec X,Vec F,void * ptr)124 PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
125 {
126 DM da;
127 PetscInt i, j, Mx, My, xs, ys, xm, ym;
128 PetscReal two = 2.0, hx, hy, sx, sy;
129 PetscScalar u, uxx, uyy, **x, **f;
130 Vec localX;
131
132 PetscFunctionBeginUser;
133 PetscCall(TSGetDM(ts, &da));
134 PetscCall(DMGetLocalVector(da, &localX));
135 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));
136
137 hx = 1.0 / (PetscReal)(Mx - 1);
138 sx = 1.0 / (hx * hx);
139 hy = 1.0 / (PetscReal)(My - 1);
140 sy = 1.0 / (hy * hy);
141
142 /*
143 Scatter ghost points to local vector,using the 2-step process
144 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
145 By placing code between these two statements, computations can be
146 done while messages are in transition.
147 */
148 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
149 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));
150
151 /*
152 Get pointers to vector data
153 */
154 PetscCall(DMDAVecGetArrayRead(da, localX, &x));
155 PetscCall(DMDAVecGetArray(da, F, &f));
156
157 /*
158 Get local grid boundaries
159 */
160 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
161
162 /*
163 Compute function over the locally owned part of the grid
164 */
165 for (j = ys; j < ys + ym; j++) {
166 for (i = xs; i < xs + xm; i++) {
167 if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
168 f[j][i] = x[j][i];
169 continue;
170 }
171 u = x[j][i];
172 uxx = (two * u - x[j][i - 1] - x[j][i + 1]) * sx;
173 uyy = (two * u - x[j - 1][i] - x[j + 1][i]) * sy;
174 /* f[j][i] = -(uxx + uyy); */
175 f[j][i] = -u * (uxx + uyy) - (4.0 - 1.0) * ((x[j][i + 1] - x[j][i - 1]) * (x[j][i + 1] - x[j][i - 1]) * .25 * sx + (x[j + 1][i] - x[j - 1][i]) * (x[j + 1][i] - x[j - 1][i]) * .25 * sy);
176 }
177 }
178
179 /*
180 Restore vectors
181 */
182 PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
183 PetscCall(DMDAVecRestoreArray(da, F, &f));
184 PetscCall(DMRestoreLocalVector(da, &localX));
185 PetscCall(PetscLogFlops(11.0 * ym * xm));
186 PetscFunctionReturn(PETSC_SUCCESS);
187 }
188
189 /* ------------------------------------------------------------------- */
FormInitialSolution(DM da,Vec U)190 PetscErrorCode FormInitialSolution(DM da, Vec U)
191 {
192 PetscInt i, j, xs, ys, xm, ym, Mx, My;
193 PetscScalar **u;
194 PetscReal hx, hy, x, y, r;
195
196 PetscFunctionBeginUser;
197 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));
198
199 hx = 1.0 / (PetscReal)(Mx - 1);
200 hy = 1.0 / (PetscReal)(My - 1);
201
202 /*
203 Get pointers to vector data
204 */
205 PetscCall(DMDAVecGetArray(da, U, &u));
206
207 /*
208 Get local grid boundaries
209 */
210 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
211
212 /*
213 Compute function over the locally owned part of the grid
214 */
215 for (j = ys; j < ys + ym; j++) {
216 y = j * hy;
217 for (i = xs; i < xs + xm; i++) {
218 x = i * hx;
219 r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
220 if (r < .125) u[j][i] = PetscExpReal(-30.0 * r * r * r);
221 else u[j][i] = 0.0;
222 }
223 }
224
225 /*
226 Restore vectors
227 */
228 PetscCall(DMDAVecRestoreArray(da, U, &u));
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
MyTSMonitor(TS ts,PetscInt step,PetscReal ptime,Vec v,PetscCtx ctx)232 PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, PetscCtx ctx)
233 {
234 PetscReal norm;
235 MPI_Comm comm;
236
237 PetscFunctionBeginUser;
238 if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* step of -1 indicates an interpolated solution */
239 PetscCall(VecNorm(v, NORM_2, &norm));
240 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
241 PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
242 PetscFunctionReturn(PETSC_SUCCESS);
243 }
244
245 /*
246 MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
247 Input Parameters:
248 snes - the SNES context
249 its - iteration number
250 fnorm - 2-norm function value (may be estimated)
251 ctx - optional user-defined context for private data for the
252 monitor routine, as set by SNESMonitorSet()
253 */
MySNESMonitor(SNES snes,PetscInt its,PetscReal fnorm,PetscViewerAndFormat * vf)254 PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
255 {
256 PetscFunctionBeginUser;
257 PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf));
258 PetscFunctionReturn(PETSC_SUCCESS);
259 }
260
261 /*TEST
262
263 test:
264 args: -ts_max_steps 5
265
266 test:
267 suffix: 2
268 args: -ts_max_steps 5 -snes_mf_operator
269
270 test:
271 suffix: 3
272 args: -ts_max_steps 5 -snes_mf -pc_type none
273
274 TEST*/
275