Actual source code: ex17.c
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
8: Boundary conditions:
9: Dirichlet BC:
10: At x=0, x=1, u = 0.0
12: Neumann BC:
13: At x=0, x=1: du(x,t)/dx = 0
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: */
22: #include <petscdm.h>
23: #include <petscdmda.h>
24: #include <petscts.h>
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};
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;
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 *);
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;
58: PetscFunctionBeginUser;
59: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
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));
68: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
69: Extract global vectors from DMDA;
70: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
71: PetscCall(DMCreateGlobalVector(da, &u));
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;
78: PetscCall(PetscOptionsGetInt(NULL, NULL, "-boundary", &user.boundary, NULL));
79: PetscCall(PetscOptionsHasName(NULL, NULL, "-viewJacobian", &user.viewJacobian));
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));
90: PetscCall(DMSetMatType(da, MATAIJ));
91: PetscCall(DMCreateMatrix(da, &J));
92: jacType = JACOBIAN_ANALYTIC; /* use user-provide Jacobian */
94: PetscCall(TSSetDM(ts, da)); /* Use TSGetDM() to access. Setting here allows easy use of geometric multigrid. */
96: ftime = 1.0;
97: PetscCall(TSSetMaxTime(ts, ftime));
98: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
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));
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: }
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Set runtime options
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: PetscCall(TSSetFromOptions(ts));
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Integrate ODE system
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: PetscCall(TSSolve(ts, u));
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));
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: /* ------------------------------------------------------------------- */
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;
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));
170: hx = 1.0 / (PetscReal)(Mx - 1);
171: sx = 1.0 / (hx * hx);
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));
182: /* Get pointers to vector data */
183: PetscCall(DMDAVecGetArrayRead(da, localU, &u));
184: PetscCall(DMDAVecGetArrayRead(da, Udot, &udot));
185: PetscCall(DMDAVecGetArray(da, F, &f));
187: /* Get local grid boundaries */
188: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
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: }
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: }
210: /* --------------------------------------------------------------------- */
211: /*
212: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
213: */
214: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat J, Mat Jpre, void *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];
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: }
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: }
270: /* ------------------------------------------------------------------- */
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;
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));
284: hx = 1.0 / (PetscReal)(Mx - 1);
286: /* Get pointers to vector data */
287: PetscCall(DMDAVecGetArray(da, U, &u));
289: /* Get local grid boundaries */
290: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
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: }
300: /* Restore vectors */
301: PetscCall(DMDAVecRestoreArray(da, U, &u));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: /*TEST
307: test:
308: requires: !single
309: args: -da_grid_x 40 -ts_max_steps 2 -snes_monitor_short -ksp_monitor_short -ts_monitor
311: test:
312: suffix: 2
313: requires: !single
314: args: -da_grid_x 100 -ts_type theta -ts_theta_theta 0.5
316: TEST*/