1 static char help[] = "Solves a time-dependent linear PDE with discontinuous right-hand side.\n";
2
3 /* ------------------------------------------------------------------------
4
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)
20
21 Reference: L. Shampine and S. Thompson, "Event Location for Ordinary Differential Equations",
22 http://www.radford.edu/~thompson/webddes/eventsweb.pdf
23 ------------------------------------------------------------------------- */
24
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;
36
InitialConditions(Vec U,DM da,AppCtx * app)37 PetscErrorCode InitialConditions(Vec U, DM da, AppCtx *app)
38 {
39 Vec xcoord;
40 PetscScalar *x, *u;
41 PetscInt lsize, M, xs, xm, i;
42
43 PetscFunctionBeginUser;
44 PetscCall(DMGetCoordinates(da, &xcoord));
45 PetscCall(DMDAVecGetArrayRead(da, xcoord, &x));
46
47 PetscCall(VecGetLocalSize(U, &lsize));
48 PetscCall(PetscMalloc1(lsize, &app->sw));
49
50 PetscCall(DMDAVecGetArray(da, U, &u));
51
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));
54
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;
59
60 app->sw[i - xs] = 1;
61 }
62 PetscCall(DMDAVecRestoreArray(da, U, &u));
63 PetscCall(DMDAVecRestoreArrayRead(da, xcoord, &x));
64 PetscFunctionReturn(PETSC_SUCCESS);
65 }
66
EventFunction(TS ts,PetscReal t,Vec U,PetscReal * fvalue,PetscCtx ctx)67 PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, PetscCtx ctx)
68 {
69 AppCtx *app = (AppCtx *)ctx;
70 const PetscScalar *u;
71 PetscInt i, lsize;
72
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 }
80
PostEventFunction(TS ts,PetscInt nevents_zero,PetscInt events_zero[],PetscReal t,Vec U,PetscBool forwardsolve,PetscCtx ctx)81 PetscErrorCode PostEventFunction(TS ts, PetscInt nevents_zero, PetscInt events_zero[], PetscReal t, Vec U, PetscBool forwardsolve, PetscCtx ctx)
82 {
83 AppCtx *app = (AppCtx *)ctx;
84 PetscInt i, idx;
85
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 }
93
94 /*
95 Defines the ODE passed to the ODE solver
96 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)97 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx 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;
106
107 PetscFunctionBeginUser;
108 PetscCall(TSGetDM(ts, &da));
109
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));
112
113 PetscCall(DMGetLocalVector(da, &Ulocal));
114 PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, Ulocal));
115 PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, Ulocal));
116
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));
122
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 }
132
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 }
139
140 /*
141 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
142 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,PetscCtx ctx)143 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx 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;
151
152 PetscFunctionBeginUser;
153 PetscCall(TSGetDM(ts, &da));
154
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));
157
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 }
185
main(int argc,char ** argv)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;
194
195 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196 Initialize program
197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198 PetscFunctionBeginUser;
199 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
200
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();
209
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));
219
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));
229
230 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231 Set initial conditions
232 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
233 PetscCall(TSSetSolution(ts, U));
234
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));
243
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));
256
257 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258 Run timestepping solver
259 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260 PetscCall(TSSolve(ts, U));
261
262 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
263 Free work space. All PETSc objects should be destroyed when they are no longer needed.
264 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
265
266 PetscCall(MatDestroy(&J));
267 PetscCall(VecDestroy(&U));
268 PetscCall(DMDestroy(&da));
269 PetscCall(TSDestroy(&ts));
270 PetscCall(PetscFree(direction));
271 PetscCall(PetscFree(terminate));
272
273 PetscCall(PetscFree(app.sw));
274 PetscCall(PetscFinalize());
275 return 0;
276 }
277