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;
45: DMGetCoordinates(da, &xcoord);
46: DMDAVecGetArrayRead(da, xcoord, &x);
48: VecGetLocalSize(U, &lsize);
49: PetscMalloc1(lsize, &app->sw);
51: DMDAVecGetArray(da, U, &u);
53: DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
54: 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: DMDAVecRestoreArray(da, U, &u);
64: DMDAVecRestoreArrayRead(da, xcoord, &x);
65: return 0;
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;
75: VecGetLocalSize(U, &lsize);
76: VecGetArrayRead(U, &u);
77: for (i = 0; i < lsize; i++) fvalue[i] = u[i] - app->uc;
78: VecRestoreArrayRead(U, &u);
79: return 0;
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;
88: for (i = 0; i < nevents_zero; i++) {
89: idx = events_zero[i];
90: app->sw[idx] = 0;
91: }
92: return 0;
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;
109: TSGetDM(ts, &da);
111: DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
112: DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
114: DMGetLocalVector(da, &Ulocal);
115: DMGlobalToLocalBegin(da, U, INSERT_VALUES, Ulocal);
116: DMGlobalToLocalEnd(da, U, INSERT_VALUES, Ulocal);
118: h = 1.0 / (M - 1);
119: h2 = h * h;
120: DMDAVecGetArrayRead(da, Udot, &udot);
121: DMDAVecGetArrayRead(da, Ulocal, &u);
122: 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: DMDAVecRestoreArrayRead(da, Udot, &udot);
135: DMDAVecRestoreArrayRead(da, Ulocal, &u);
136: DMDAVecRestoreArray(da, F, &f);
137: DMRestoreLocalVector(da, &Ulocal);
139: return 0;
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;
155: TSGetDM(ts, &da);
157: DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
158: 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: 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: 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: MatSetValuesStencil(A, 1, &row, 3, col, v, INSERT_VALUES);
181: }
182: }
183: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
184: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
185: return 0;
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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: PetscInitialize(&argc, &argv, (char *)0, help);
203: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex22 options", "");
204: {
205: app.A = 200000;
206: PetscOptionsReal("-A", "", "", app.A, &app.A, NULL);
207: app.uc = 0.5;
208: PetscOptionsReal("-uc", "", "", app.uc, &app.uc, NULL);
209: }
210: PetscOptionsEnd();
212: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, -n, 1, 1, 0, &da);
213: DMSetFromOptions(da);
214: DMSetUp(da);
215: DMDASetUniformCoordinates(da, 0.0, 1.0, 0, 0, 0, 0);
216: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: Create necessary matrix and vectors
218: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
219: DMCreateMatrix(da, &J);
220: DMCreateGlobalVector(da, &U);
222: InitialConditions(U, da, &app);
223: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
224: Create timestepping solver context
225: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
226: TSCreate(PETSC_COMM_WORLD, &ts);
227: TSSetProblemType(ts, TS_NONLINEAR);
228: TSSetType(ts, TSROSW);
229: TSSetIFunction(ts, NULL, (TSIFunction)IFunction, (void *)&app);
230: TSSetIJacobian(ts, J, J, (TSIJacobian)IJacobian, (void *)&app);
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Set initial conditions
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: TSSetSolution(ts, U);
237: TSSetDM(ts, da);
238: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: Set solver options
240: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
241: TSSetTimeStep(ts, 0.1);
242: TSSetMaxTime(ts, 30.0);
243: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
244: TSSetFromOptions(ts);
246: PetscInt lsize;
247: VecGetLocalSize(U, &lsize);
248: PetscInt *direction;
249: PetscBool *terminate;
250: PetscInt i;
251: PetscMalloc1(lsize, &direction);
252: PetscMalloc1(lsize, &terminate);
253: for (i = 0; i < lsize; i++) {
254: direction[i] = -1;
255: terminate[i] = PETSC_FALSE;
256: }
257: TSSetEventHandler(ts, lsize, direction, terminate, EventFunction, PostEventFunction, (void *)&app);
259: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: Run timestepping solver
261: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
262: TSSolve(ts, U);
264: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
265: Free work space. All PETSc objects should be destroyed when they are no longer needed.
266: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
268: MatDestroy(&J);
269: VecDestroy(&U);
270: DMDestroy(&da);
271: TSDestroy(&ts);
272: PetscFree(direction);
273: PetscFree(terminate);
275: PetscFree(app.sw);
276: PetscFinalize();
277: return 0;
278: }