1 static const char help[] = "Time-dependent PDE in 1d. Simplified from ex15.c for illustrating how to solve DAEs. \n";
2 /*
3 u_t = uxx
4 0 < x < 1;
5 At t=0: u(x) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5)) < .125
6 u(x) = 0.0 if r >= .125
7
8 Boundary conditions:
9 Dirichlet BC:
10 At x=0, x=1, u = 0.0
11
12 Neumann BC:
13 At x=0, x=1: du(x,t)/dx = 0
14
15 mpiexec -n 2 ./ex17 -da_grid_x 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
16 ./ex17 -da_grid_x 40 -monitor_solution
17 ./ex17 -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5 # Midpoint is not L-stable
18 ./ex17 -jac_type fd_coloring -da_grid_x 500 -boundary 1
19 ./ex17 -da_grid_x 100 -ts_type gl -ts_adapt_type none -ts_max_steps 2
20 */
21
22 #include <petscdm.h>
23 #include <petscdmda.h>
24 #include <petscts.h>
25
26 typedef enum {
27 JACOBIAN_ANALYTIC,
28 JACOBIAN_FD_COLORING,
29 JACOBIAN_FD_FULL
30 } JacobianType;
31 static const char *const JacobianTypes[] = {"analytic", "fd_coloring", "fd_full", "JacobianType", "fd_", 0};
32
33 /*
34 User-defined data structures and routines
35 */
36 typedef struct {
37 PetscReal c;
38 PetscInt boundary; /* Type of boundary condition */
39 PetscBool viewJacobian;
40 } AppCtx;
41
42 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
43 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
44 static PetscErrorCode FormInitialSolution(TS, Vec, void *);
45
main(int argc,char ** argv)46 int main(int argc, char **argv)
47 {
48 TS ts; /* nonlinear solver */
49 Vec u; /* solution, residual vectors */
50 Mat J; /* Jacobian matrix */
51 PetscInt nsteps;
52 PetscReal vmin, vmax, norm;
53 DM da;
54 PetscReal ftime, dt;
55 AppCtx user; /* user-defined work context */
56 JacobianType jacType;
57
58 PetscFunctionBeginUser;
59 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
60
61 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62 Create distributed array (DMDA) to manage parallel grid and vectors
63 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 1, 1, NULL, &da));
65 PetscCall(DMSetFromOptions(da));
66 PetscCall(DMSetUp(da));
67
68 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69 Extract global vectors from DMDA;
70 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71 PetscCall(DMCreateGlobalVector(da, &u));
72
73 /* Initialize user application context */
74 user.c = -30.0;
75 user.boundary = 0; /* 0: Dirichlet BC; 1: Neumann BC */
76 user.viewJacobian = PETSC_FALSE;
77
78 PetscCall(PetscOptionsGetInt(NULL, NULL, "-boundary", &user.boundary, NULL));
79 PetscCall(PetscOptionsHasName(NULL, NULL, "-viewJacobian", &user.viewJacobian));
80
81 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 Create timestepping solver context
83 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
85 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
86 PetscCall(TSSetType(ts, TSTHETA));
87 PetscCall(TSThetaSetTheta(ts, 1.0)); /* Make the Theta method behave like backward Euler */
88 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
89
90 PetscCall(DMSetMatType(da, MATAIJ));
91 PetscCall(DMCreateMatrix(da, &J));
92 jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */
93
94 PetscCall(TSSetDM(ts, da)); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */
95
96 ftime = 1.0;
97 PetscCall(TSSetMaxTime(ts, ftime));
98 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
99
100 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101 Set initial conditions
102 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103 PetscCall(FormInitialSolution(ts, u, &user));
104 PetscCall(TSSetSolution(ts, u));
105 dt = .01;
106 PetscCall(TSSetTimeStep(ts, dt));
107
108 /* Use slow fd Jacobian or fast fd Jacobian with colorings.
109 Note: this requires snes which is not created until TSSetUp()/TSSetFromOptions() is called */
110 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Options for Jacobian evaluation", NULL);
111 PetscCall(PetscOptionsEnum("-jac_type", "Type of Jacobian", "", JacobianTypes, (PetscEnum)jacType, (PetscEnum *)&jacType, 0));
112 PetscOptionsEnd();
113 if (jacType == JACOBIAN_ANALYTIC) {
114 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));
115 } else if (jacType == JACOBIAN_FD_COLORING) {
116 SNES snes;
117 PetscCall(TSGetSNES(ts, &snes));
118 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, 0));
119 } else if (jacType == JACOBIAN_FD_FULL) {
120 SNES snes;
121 PetscCall(TSGetSNES(ts, &snes));
122 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, &user));
123 }
124
125 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 Set runtime options
127 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128 PetscCall(TSSetFromOptions(ts));
129
130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131 Integrate ODE system
132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133 PetscCall(TSSolve(ts, u));
134
135 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136 Compute diagnostics of the solution
137 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138 PetscCall(VecNorm(u, NORM_1, &norm));
139 PetscCall(VecMax(u, NULL, &vmax));
140 PetscCall(VecMin(u, NULL, &vmin));
141 PetscCall(TSGetStepNumber(ts, &nsteps));
142 PetscCall(TSGetTime(ts, &ftime));
143 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "timestep %" PetscInt_FMT ": time %g, solution norm %g, max %g, min %g\n", nsteps, (double)ftime, (double)norm, (double)vmax, (double)vmin));
144
145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146 Free work space.
147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148 PetscCall(MatDestroy(&J));
149 PetscCall(VecDestroy(&u));
150 PetscCall(TSDestroy(&ts));
151 PetscCall(DMDestroy(&da));
152 PetscCall(PetscFinalize());
153 return 0;
154 }
155 /* ------------------------------------------------------------------- */
FormIFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)156 static PetscErrorCode FormIFunction(TS ts, PetscReal ftime, Vec U, Vec Udot, Vec F, void *ptr)
157 {
158 AppCtx *user = (AppCtx *)ptr;
159 DM da;
160 PetscInt i, Mx, xs, xm;
161 PetscReal hx, sx;
162 PetscScalar *u, *udot, *f;
163 Vec localU;
164
165 PetscFunctionBeginUser;
166 PetscCall(TSGetDM(ts, &da));
167 PetscCall(DMGetLocalVector(da, &localU));
168 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
169
170 hx = 1.0 / (PetscReal)(Mx - 1);
171 sx = 1.0 / (hx * hx);
172
173 /*
174 Scatter ghost points to local vector,using the 2-step process
175 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
176 By placing code between these two statements, computations can be
177 done while messages are in transition.
178 */
179 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
180 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));
181
182 /* Get pointers to vector data */
183 PetscCall(DMDAVecGetArrayRead(da, localU, &u));
184 PetscCall(DMDAVecGetArrayRead(da, Udot, &udot));
185 PetscCall(DMDAVecGetArray(da, F, &f));
186
187 /* Get local grid boundaries */
188 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
189
190 /* Compute function over the locally owned part of the grid */
191 for (i = xs; i < xs + xm; i++) {
192 if (user->boundary == 0) { /* Dirichlet BC */
193 if (i == 0 || i == Mx - 1) f[i] = u[i]; /* F = U */
194 else f[i] = udot[i] + (2. * u[i] - u[i - 1] - u[i + 1]) * sx;
195 } else { /* Neumann BC */
196 if (i == 0) f[i] = u[0] - u[1];
197 else if (i == Mx - 1) f[i] = u[i] - u[i - 1];
198 else f[i] = udot[i] + (2. * u[i] - u[i - 1] - u[i + 1]) * sx;
199 }
200 }
201
202 /* Restore vectors */
203 PetscCall(DMDAVecRestoreArrayRead(da, localU, &u));
204 PetscCall(DMDAVecRestoreArrayRead(da, Udot, &udot));
205 PetscCall(DMDAVecRestoreArray(da, F, &f));
206 PetscCall(DMRestoreLocalVector(da, &localU));
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
210 /* --------------------------------------------------------------------- */
211 /*
212 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
213 */
FormIJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat J,Mat Jpre,PetscCtx ctx)214 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, PetscCtx ctx)
215 {
216 PetscInt i, rstart, rend, Mx;
217 PetscReal hx, sx;
218 AppCtx *user = (AppCtx *)ctx;
219 DM da;
220 MatStencil col[3], row;
221 PetscInt nc;
222 PetscScalar vals[3];
223
224 PetscFunctionBeginUser;
225 PetscCall(TSGetDM(ts, &da));
226 PetscCall(MatGetOwnershipRange(Jpre, &rstart, &rend));
227 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
228 hx = 1.0 / (PetscReal)(Mx - 1);
229 sx = 1.0 / (hx * hx);
230 for (i = rstart; i < rend; i++) {
231 nc = 0;
232 row.i = i;
233 if (user->boundary == 0 && (i == 0 || i == Mx - 1)) {
234 col[nc].i = i;
235 vals[nc++] = 1.0;
236 } else if (user->boundary > 0 && i == 0) { /* Left Neumann */
237 col[nc].i = i;
238 vals[nc++] = 1.0;
239 col[nc].i = i + 1;
240 vals[nc++] = -1.0;
241 } else if (user->boundary > 0 && i == Mx - 1) { /* Right Neumann */
242 col[nc].i = i - 1;
243 vals[nc++] = -1.0;
244 col[nc].i = i;
245 vals[nc++] = 1.0;
246 } else { /* Interior */
247 col[nc].i = i - 1;
248 vals[nc++] = -1.0 * sx;
249 col[nc].i = i;
250 vals[nc++] = 2.0 * sx + a;
251 col[nc].i = i + 1;
252 vals[nc++] = -1.0 * sx;
253 }
254 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, vals, INSERT_VALUES));
255 }
256
257 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
258 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
259 if (J != Jpre) {
260 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
261 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
262 }
263 if (user->viewJacobian) {
264 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Jpre:\n"));
265 PetscCall(MatView(Jpre, PETSC_VIEWER_STDOUT_WORLD));
266 }
267 PetscFunctionReturn(PETSC_SUCCESS);
268 }
269
270 /* ------------------------------------------------------------------- */
FormInitialSolution(TS ts,Vec U,void * ptr)271 PetscErrorCode FormInitialSolution(TS ts, Vec U, void *ptr)
272 {
273 AppCtx *user = (AppCtx *)ptr;
274 PetscReal c = user->c;
275 DM da;
276 PetscInt i, xs, xm, Mx;
277 PetscScalar *u;
278 PetscReal hx, x, r;
279
280 PetscFunctionBeginUser;
281 PetscCall(TSGetDM(ts, &da));
282 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE));
283
284 hx = 1.0 / (PetscReal)(Mx - 1);
285
286 /* Get pointers to vector data */
287 PetscCall(DMDAVecGetArray(da, U, &u));
288
289 /* Get local grid boundaries */
290 PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
291
292 /* Compute function over the locally owned part of the grid */
293 for (i = xs; i < xs + xm; i++) {
294 x = i * hx;
295 r = PetscSqrtReal((x - .5) * (x - .5));
296 if (r < .125) u[i] = PetscExpReal(c * r * r * r);
297 else u[i] = 0.0;
298 }
299
300 /* Restore vectors */
301 PetscCall(DMDAVecRestoreArray(da, U, &u));
302 PetscFunctionReturn(PETSC_SUCCESS);
303 }
304
305 /*TEST
306
307 test:
308 requires: !single
309 args: -da_grid_x 40 -ts_max_steps 2 -snes_monitor_short -ksp_monitor_short -ts_monitor
310
311 test:
312 suffix: 2
313 requires: !single
314 args: -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5
315
316 TEST*/
317