Actual source code: ex32.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, PetscReal *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] = PetscRealPart(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));
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: /*
    141:      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
    142: */
    143: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *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;

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

    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));

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

    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;

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

    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();

    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));

    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));

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

    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));

    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));

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

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

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

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