Actual source code: ex22.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, PetscScalar *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] = 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));

138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: /*
142:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
143: */
144: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
145: {
146:   AppCtx     *app = (AppCtx *)ctx;
147:   DM          da;
148:   MatStencil  row, col[3];
149:   PetscScalar v[3];
150:   PetscInt    M, xs, xm, i;
151:   PetscReal   h, h2;

153:   PetscFunctionBeginUser;
154:   PetscCall(TSGetDM(ts, &da));

156:   PetscCall(DMDAGetInfo(da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
157:   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));

159:   h  = 1.0 / (M - 1);
160:   h2 = h * h;
161:   for (i = xs; i < xs + xm; i++) {
162:     row.i = i;
163:     if (i == 0) {
164:       v[0] = 1.0;
165:       PetscCall(MatSetValuesStencil(A, 1, &row, 1, &row, v, INSERT_VALUES));
166:     } else if (i == M - 1) {
167:       col[0].i = i;
168:       v[0]     = 1 / h;
169:       col[1].i = i - 1;
170:       v[1]     = -1 / h;
171:       PetscCall(MatSetValuesStencil(A, 1, &row, 2, col, v, INSERT_VALUES));
172:     } else {
173:       col[0].i = i + 1;
174:       v[0]     = 1 / h2;
175:       col[1].i = i;
176:       v[1]     = -2 / h2 + app->sw[i - xs] * (-app->A) - a;
177:       col[2].i = i - 1;
178:       v[2]     = 1 / h2;
179:       PetscCall(MatSetValuesStencil(A, 1, &row, 3, col, v, INSERT_VALUES));
180:     }
181:   }
182:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
183:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: int main(int argc, char **argv)
188: {
189:   TS       ts; /* ODE integrator */
190:   Vec      U;  /* solution will be stored here */
191:   Mat      J;  /* Jacobian matrix */
192:   PetscInt n = 16;
193:   AppCtx   app;
194:   DM       da;

196:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197:      Initialize program
198:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199:   PetscFunctionBeginUser;
200:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

202:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex22 options", "");
203:   {
204:     app.A = 200000;
205:     PetscCall(PetscOptionsReal("-A", "", "", app.A, &app.A, NULL));
206:     app.uc = 0.5;
207:     PetscCall(PetscOptionsReal("-uc", "", "", app.uc, &app.uc, NULL));
208:   }
209:   PetscOptionsEnd();

211:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, -n, 1, 1, 0, &da));
212:   PetscCall(DMSetFromOptions(da));
213:   PetscCall(DMSetUp(da));
214:   PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0, 0, 0, 0));
215:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216:     Create necessary matrix and vectors
217:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218:   PetscCall(DMCreateMatrix(da, &J));
219:   PetscCall(DMCreateGlobalVector(da, &U));

221:   PetscCall(InitialConditions(U, da, &app));
222:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
223:      Create timestepping solver context
224:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
226:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
227:   PetscCall(TSSetType(ts, TSROSW));
228:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, (void *)&app));
229:   PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobian)IJacobian, (void *)&app));

231:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232:      Set initial conditions
233:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234:   PetscCall(TSSetSolution(ts, U));

236:   PetscCall(TSSetDM(ts, da));
237:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238:      Set solver options
239:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
240:   PetscCall(TSSetTimeStep(ts, 0.1));
241:   PetscCall(TSSetMaxTime(ts, 30.0));
242:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
243:   PetscCall(TSSetFromOptions(ts));

245:   PetscInt lsize;
246:   PetscCall(VecGetLocalSize(U, &lsize));
247:   PetscInt  *direction;
248:   PetscBool *terminate;
249:   PetscInt   i;
250:   PetscCall(PetscMalloc1(lsize, &direction));
251:   PetscCall(PetscMalloc1(lsize, &terminate));
252:   for (i = 0; i < lsize; i++) {
253:     direction[i] = -1;
254:     terminate[i] = PETSC_FALSE;
255:   }
256:   PetscCall(TSSetEventHandler(ts, lsize, direction, terminate, EventFunction, PostEventFunction, (void *)&app));

258:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259:      Run timestepping solver
260:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
261:   PetscCall(TSSolve(ts, U));

263:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
264:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
265:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

267:   PetscCall(MatDestroy(&J));
268:   PetscCall(VecDestroy(&U));
269:   PetscCall(DMDestroy(&da));
270:   PetscCall(TSDestroy(&ts));
271:   PetscCall(PetscFree(direction));
272:   PetscCall(PetscFree(terminate));

274:   PetscCall(PetscFree(app.sw));
275:   PetscCall(PetscFinalize());
276:   return 0;
277: }