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