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: }