Actual source code: ex25.c

      1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d formulated as a PDAE. Demonstrates solving PDEs with algebraic constraints (PDAE).\n";
      2: /*
      3:    u_t - alpha u_xx = A + u^2 v - (B+1) u
      4:    v_t - alpha v_xx = B u - u^2 v
      5:    0 < x < 1;
      6:    A = 1, B = 3, alpha = 1/50

      8:    Initial conditions:
      9:    u(x,0) = 1 + sin(2 pi x)
     10:    v(x,0) = 3

     12:    Boundary conditions:
     13:    u(0,t) = u(1,t) = 1
     14:    v(0,t) = v(1,t) = 3
     15: */

     17: #include <petscdm.h>
     18: #include <petscdmda.h>
     19: #include <petscts.h>

     21: typedef struct {
     22:   PetscScalar u, v;
     23: } Field;

     25: typedef struct _User *User;
     26: struct _User {
     27:   PetscReal A, B;          /* Reaction coefficients */
     28:   PetscReal alpha;         /* Diffusion coefficient */
     29:   PetscReal uleft, uright; /* Dirichlet boundary conditions */
     30:   PetscReal vleft, vright; /* Dirichlet boundary conditions */
     31: };

     33: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
     34: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     35: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
     36: static PetscErrorCode FormInitialSolution(TS, Vec, void *);

     38: int main(int argc, char **argv)
     39: {
     40:   TS                ts; /* nonlinear solver */
     41:   Vec               X;  /* solution, residual vectors */
     42:   Mat               J;  /* Jacobian matrix */
     43:   PetscInt          steps, mx;
     44:   DM                da;
     45:   PetscReal         ftime, hx, dt;
     46:   struct _User      user; /* user-defined work context */
     47:   TSConvergedReason reason;

     49:   PetscFunctionBeginUser;
     50:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     52:      Create distributed array (DMDA) to manage parallel grid and vectors
     53:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     54:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da));
     55:   PetscCall(DMSetFromOptions(da));
     56:   PetscCall(DMSetUp(da));

     58:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     59:      Extract global vectors from DMDA;
     60:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     61:   PetscCall(DMCreateGlobalVector(da, &X));

     63:   /* Initialize user application context */
     64:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
     65:   {
     66:     user.A      = 1;
     67:     user.B      = 3;
     68:     user.alpha  = 0.02;
     69:     user.uleft  = 1;
     70:     user.uright = 1;
     71:     user.vleft  = 3;
     72:     user.vright = 3;
     73:     PetscCall(PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL));
     74:     PetscCall(PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL));
     75:     PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL));
     76:     PetscCall(PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL));
     77:     PetscCall(PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL));
     78:     PetscCall(PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL));
     79:     PetscCall(PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL));
     80:   }
     81:   PetscOptionsEnd();

     83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     84:      Create timestepping solver context
     85:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     86:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     87:   PetscCall(TSSetDM(ts, da));
     88:   PetscCall(TSSetType(ts, TSARKIMEX));
     89:   PetscCall(TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1));
     90:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
     91:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
     92:   PetscCall(DMSetMatType(da, MATAIJ));
     93:   PetscCall(DMCreateMatrix(da, &J));
     94:   PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));

     96:   ftime = 10.0;
     97:   PetscCall(TSSetMaxTime(ts, ftime));
     98:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101:      Set initial conditions
    102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    103:   PetscCall(FormInitialSolution(ts, X, &user));
    104:   PetscCall(TSSetSolution(ts, X));
    105:   PetscCall(VecGetSize(X, &mx));
    106:   hx = 1.0 / (PetscReal)(mx / 2 - 1);
    107:   dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
    108:   PetscCall(TSSetTimeStep(ts, dt));

    110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    111:      Set runtime options
    112:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    113:   PetscCall(TSSetFromOptions(ts));

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:      Solve nonlinear system
    117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    118:   PetscCall(TSSolve(ts, X));
    119:   PetscCall(TSGetSolveTime(ts, &ftime));
    120:   PetscCall(TSGetStepNumber(ts, &steps));
    121:   PetscCall(TSGetConvergedReason(ts, &reason));
    122:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

    124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    125:      Free work space.
    126:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    127:   PetscCall(MatDestroy(&J));
    128:   PetscCall(VecDestroy(&X));
    129:   PetscCall(TSDestroy(&ts));
    130:   PetscCall(DMDestroy(&da));
    131:   PetscCall(PetscFinalize());
    132:   return 0;
    133: }

    135: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
    136: {
    137:   User          user = (User)ptr;
    138:   DM            da;
    139:   DMDALocalInfo info;
    140:   PetscInt      i;
    141:   Field        *x, *xdot, *f;
    142:   PetscReal     hx;
    143:   Vec           Xloc;

    145:   PetscFunctionBeginUser;
    146:   PetscCall(TSGetDM(ts, &da));
    147:   PetscCall(DMDAGetLocalInfo(da, &info));
    148:   hx = 1.0 / (PetscReal)(info.mx - 1);

    150:   /*
    151:      Scatter ghost points to local vector,using the 2-step process
    152:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    153:      By placing code between these two statements, computations can be
    154:      done while messages are in transition.
    155:   */
    156:   PetscCall(DMGetLocalVector(da, &Xloc));
    157:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
    158:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));

    160:   /* Get pointers to vector data */
    161:   PetscCall(DMDAVecGetArrayRead(da, Xloc, &x));
    162:   PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
    163:   PetscCall(DMDAVecGetArray(da, F, &f));

    165:   /* Compute function over the locally owned part of the grid */
    166:   for (i = info.xs; i < info.xs + info.xm; i++) {
    167:     if (i == 0) {
    168:       f[i].u = hx * (x[i].u - user->uleft);
    169:       f[i].v = hx * (x[i].v - user->vleft);
    170:     } else if (i == info.mx - 1) {
    171:       f[i].u = hx * (x[i].u - user->uright);
    172:       f[i].v = hx * (x[i].v - user->vright);
    173:     } else {
    174:       f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
    175:       f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
    176:     }
    177:   }

    179:   /* Restore vectors */
    180:   PetscCall(DMDAVecRestoreArrayRead(da, Xloc, &x));
    181:   PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
    182:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    183:   PetscCall(DMRestoreLocalVector(da, &Xloc));
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    188: {
    189:   User          user = (User)ptr;
    190:   DM            da;
    191:   DMDALocalInfo info;
    192:   PetscInt      i;
    193:   PetscReal     hx;
    194:   Field        *x, *f;

    196:   PetscFunctionBeginUser;
    197:   PetscCall(TSGetDM(ts, &da));
    198:   PetscCall(DMDAGetLocalInfo(da, &info));
    199:   hx = 1.0 / (PetscReal)(info.mx - 1);

    201:   /* Get pointers to vector data */
    202:   PetscCall(DMDAVecGetArrayRead(da, X, &x));
    203:   PetscCall(DMDAVecGetArray(da, F, &f));

    205:   /* Compute function over the locally owned part of the grid */
    206:   for (i = info.xs; i < info.xs + info.xm; i++) {
    207:     PetscScalar u = x[i].u, v = x[i].v;
    208:     f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u);
    209:     f[i].v = hx * (user->B * u - u * u * v);
    210:   }

    212:   /* Restore vectors */
    213:   PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
    214:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: /* --------------------------------------------------------------------- */
    219: /*
    220:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
    221: */
    222: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
    223: {
    224:   User          user = (User)ptr;
    225:   DMDALocalInfo info;
    226:   PetscInt      i;
    227:   PetscReal     hx;
    228:   DM            da;
    229:   Field        *x, *xdot;

    231:   PetscFunctionBeginUser;
    232:   PetscCall(TSGetDM(ts, &da));
    233:   PetscCall(DMDAGetLocalInfo(da, &info));
    234:   hx = 1.0 / (PetscReal)(info.mx - 1);

    236:   /* Get pointers to vector data */
    237:   PetscCall(DMDAVecGetArrayRead(da, X, &x));
    238:   PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));

    240:   /* Compute function over the locally owned part of the grid */
    241:   for (i = info.xs; i < info.xs + info.xm; i++) {
    242:     if (i == 0 || i == info.mx - 1) {
    243:       const PetscInt    row = i, col = i;
    244:       const PetscScalar vals[2][2] = {
    245:         {hx, 0 },
    246:         {0,  hx}
    247:       };
    248:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES));
    249:     } else {
    250:       const PetscInt    row = i, col[] = {i - 1, i, i + 1};
    251:       const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
    252:       const PetscScalar vals[2][3][2] = {
    253:         {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
    254:         {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
    255:       };
    256:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES));
    257:     }
    258:   }

    260:   /* Restore vectors */
    261:   PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
    262:   PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));

    264:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    265:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    266:   if (J != Jpre) {
    267:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    268:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    269:   }
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
    274: {
    275:   User          user = (User)ctx;
    276:   DM            da;
    277:   PetscInt      i;
    278:   DMDALocalInfo info;
    279:   Field        *x;
    280:   PetscReal     hx;

    282:   PetscFunctionBeginUser;
    283:   PetscCall(TSGetDM(ts, &da));
    284:   PetscCall(DMDAGetLocalInfo(da, &info));
    285:   hx = 1.0 / (PetscReal)(info.mx - 1);

    287:   /* Get pointers to vector data */
    288:   PetscCall(DMDAVecGetArray(da, X, &x));

    290:   /* Compute function over the locally owned part of the grid */
    291:   for (i = info.xs; i < info.xs + info.xm; i++) {
    292:     PetscReal xi = i * hx;
    293:     x[i].u       = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi);
    294:     x[i].v       = user->vleft * (1. - xi) + user->vright * xi;
    295:   }
    296:   PetscCall(DMDAVecRestoreArray(da, X, &x));
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: /*TEST

    302:     test:
    303:       args: -nox -da_grid_x 20 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type 2p -ts_time_step 5e-2 -ts_adapt_type none

    305: TEST*/