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