Actual source code: ex32.c
1: static char help[] = "Solves a time-dependent linear PDE with discontinuous right-hand side.\n";
3: /* ------------------------------------------------------------------------
5: This program solves the one-dimensional quench front problem modeling a cooled
6: liquid rising on a hot metal rod
7: u_t = u_xx + g(u),
8: with
9: g(u) = -Au if u <= u_c,
10: = 0 if u > u_c
11: on the domain 0 <= x <= 1, with the boundary conditions
12: u(t,0) = 0, u_x(t,1) = 0,
13: and the initial condition
14: u(0,x) = 0 if 0 <= x <= 0.1,
15: = (x - 0.1)/0.15 if 0.1 < x < 0.25
16: = 1 if 0.25 <= x <= 1
17: We discretize the right-hand side using finite differences with
18: uniform grid spacing h:
19: u_xx = (u_{i+1} - 2u_{i} + u_{i-1})/(h^2)
21: Reference: L. Shampine and S. Thompson, "Event Location for Ordinary Differential Equations",
22: http://www.radford.edu/~thompson/webddes/eventsweb.pdf
23: ------------------------------------------------------------------------- */
25: #include <petscdmda.h>
26: #include <petscts.h>
27: /*
28: User-defined application context - contains data needed by the
29: application-provided call-back routines.
30: */
31: typedef struct {
32: PetscReal A;
33: PetscReal uc;
34: PetscInt *sw;
35: } AppCtx;
37: PetscErrorCode InitialConditions(Vec U, DM da, AppCtx *app)
38: {
39: Vec xcoord;
40: PetscScalar *x, *u;
41: PetscInt lsize, M, xs, xm, i;
43: PetscFunctionBeginUser;
44: PetscCall(DMGetCoordinates(da, &xcoord));
45: PetscCall(DMDAVecGetArrayRead(da, xcoord, &x));
47: PetscCall(VecGetLocalSize(U, &lsize));
48: PetscCall(PetscMalloc1(lsize, &app->sw));
50: PetscCall(DMDAVecGetArray(da, U, &u));
52: PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
53: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
55: for (i = xs; i < xs + xm; i++) {
56: if (x[i] <= 0.1) u[i] = 0.;
57: else if (x[i] > 0.1 && x[i] < 0.25) u[i] = (x[i] - 0.1) / 0.15;
58: else u[i] = 1.0;
60: app->sw[i - xs] = 1;
61: }
62: PetscCall(DMDAVecRestoreArray(da, U, &u));
63: PetscCall(DMDAVecRestoreArrayRead(da, xcoord, &x));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx)
68: {
69: AppCtx *app = (AppCtx *)ctx;
70: const PetscScalar *u;
71: PetscInt i, lsize;
73: PetscFunctionBeginUser;
74: PetscCall(VecGetLocalSize(U, &lsize));
75: PetscCall(VecGetArrayRead(U, &u));
76: for (i = 0; i < lsize; i++) fvalue[i] = PetscRealPart(u[i]) - app->uc;
77: PetscCall(VecRestoreArrayRead(U, &u));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
82: {
83: AppCtx *app = (AppCtx *)ctx;
84: PetscInt i, idx;
86: PetscFunctionBeginUser;
87: for (i = 0; i < nevents_zero; i++) {
88: idx = events_zero[i];
89: app->sw[idx] = 0;
90: }
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*
95: Defines the ODE passed to the ODE solver
96: */
97: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
98: {
99: AppCtx *app = (AppCtx *)ctx;
100: PetscScalar *f;
101: const PetscScalar *u, *udot;
102: DM da;
103: PetscInt M, xs, xm, i;
104: PetscReal h, h2;
105: Vec Ulocal;
107: PetscFunctionBeginUser;
108: PetscCall(TSGetDM(ts, &da));
110: PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
111: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
113: PetscCall(DMGetLocalVector(da, &Ulocal));
114: PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, Ulocal));
115: PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, Ulocal));
117: h = 1.0 / (M - 1);
118: h2 = h * h;
119: PetscCall(DMDAVecGetArrayRead(da, Udot, &udot));
120: PetscCall(DMDAVecGetArrayRead(da, Ulocal, &u));
121: PetscCall(DMDAVecGetArray(da, F, &f));
123: for (i = xs; i < xs + xm; i++) {
124: if (i == 0) {
125: f[i] = u[i];
126: } else if (i == M - 1) {
127: f[i] = (u[i] - u[i - 1]) / h;
128: } else {
129: f[i] = (u[i + 1] - 2 * u[i] + u[i - 1]) / h2 + app->sw[i - xs] * (-app->A * u[i]) - udot[i];
130: }
131: }
133: PetscCall(DMDAVecRestoreArrayRead(da, Udot, &udot));
134: PetscCall(DMDAVecRestoreArrayRead(da, Ulocal, &u));
135: PetscCall(DMDAVecRestoreArray(da, F, &f));
136: PetscCall(DMRestoreLocalVector(da, &Ulocal));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*
141: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
142: */
143: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
144: {
145: AppCtx *app = (AppCtx *)ctx;
146: DM da;
147: MatStencil row, col[3];
148: PetscScalar v[3];
149: PetscInt M, xs, xm, i;
150: PetscReal h, h2;
152: PetscFunctionBeginUser;
153: PetscCall(TSGetDM(ts, &da));
155: PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
156: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
158: h = 1.0 / (M - 1);
159: h2 = h * h;
160: for (i = xs; i < xs + xm; i++) {
161: row.i = i;
162: if (i == 0) {
163: v[0] = 1.0;
164: PetscCall(MatSetValuesStencil(A, 1, &row, 1, &row, v, INSERT_VALUES));
165: } else if (i == M - 1) {
166: col[0].i = i;
167: v[0] = 1 / h;
168: col[1].i = i - 1;
169: v[1] = -1 / h;
170: PetscCall(MatSetValuesStencil(A, 1, &row, 2, col, v, INSERT_VALUES));
171: } else {
172: col[0].i = i + 1;
173: v[0] = 1 / h2;
174: col[1].i = i;
175: v[1] = -2 / h2 + app->sw[i - xs] * (-app->A) - a;
176: col[2].i = i - 1;
177: v[2] = 1 / h2;
178: PetscCall(MatSetValuesStencil(A, 1, &row, 3, col, v, INSERT_VALUES));
179: }
180: }
181: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
182: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: int main(int argc, char **argv)
187: {
188: TS ts; /* ODE integrator */
189: Vec U; /* solution will be stored here */
190: Mat J; /* Jacobian matrix */
191: PetscInt n = 16;
192: AppCtx app;
193: DM da;
195: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196: Initialize program
197: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: PetscFunctionBeginUser;
199: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
201: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex22 options", "");
202: {
203: app.A = 200000;
204: PetscCall(PetscOptionsReal("-A", "", "", app.A, &app.A, NULL));
205: app.uc = 0.5;
206: PetscCall(PetscOptionsReal("-uc", "", "", app.uc, &app.uc, NULL));
207: }
208: PetscOptionsEnd();
210: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, -n, 1, 1, 0, &da));
211: PetscCall(DMSetFromOptions(da));
212: PetscCall(DMSetUp(da));
213: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0, 0, 0, 0));
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Create necessary matrix and vectors
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: PetscCall(DMCreateMatrix(da, &J));
218: PetscCall(DMCreateGlobalVector(da, &U));
220: PetscCall(InitialConditions(U, da, &app));
221: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222: Create timestepping solver context
223: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
225: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
226: PetscCall(TSSetType(ts, TSROSW));
227: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, (void *)&app));
228: PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, (void *)&app));
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Set initial conditions
232: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233: PetscCall(TSSetSolution(ts, U));
235: PetscCall(TSSetDM(ts, da));
236: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: Set solver options
238: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239: PetscCall(TSSetTimeStep(ts, 0.1));
240: PetscCall(TSSetMaxTime(ts, 30.0));
241: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
242: PetscCall(TSSetFromOptions(ts));
244: PetscInt lsize;
245: PetscCall(VecGetLocalSize(U, &lsize));
246: PetscInt *direction;
247: PetscBool *terminate;
248: PetscInt i;
249: PetscCall(PetscMalloc1(lsize, &direction));
250: PetscCall(PetscMalloc1(lsize, &terminate));
251: for (i = 0; i < lsize; i++) {
252: direction[i] = -1;
253: terminate[i] = PETSC_FALSE;
254: }
255: PetscCall(TSSetEventHandler(ts, lsize, direction, terminate, EventFunction, PostEventFunction, (void *)&app));
257: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: Run timestepping solver
259: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260: PetscCall(TSSolve(ts, U));
262: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
263: Free work space. All PETSc objects should be destroyed when they are no longer needed.
264: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
266: PetscCall(MatDestroy(&J));
267: PetscCall(VecDestroy(&U));
268: PetscCall(DMDestroy(&da));
269: PetscCall(TSDestroy(&ts));
270: PetscCall(PetscFree(direction));
271: PetscCall(PetscFree(terminate));
273: PetscCall(PetscFree(app.sw));
274: PetscCall(PetscFinalize());
275: return 0;
276: }