Actual source code: ex35.cxx

      1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and uses MOAB.\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: // PETSc includes:
     18: #include <petscts.h>
     19: #include <petscdmmoab.h>

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

     25: struct pUserCtx {
     26:   PetscReal A, B;    /* Reaction coefficients */
     27:   PetscReal alpha;   /* Diffusion coefficient */
     28:   Field     leftbc;  /* Dirichlet boundary conditions at left boundary */
     29:   Field     rightbc; /* Dirichlet boundary conditions at right boundary */
     30:   PetscInt  n, npts; /* Number of mesh points */
     31:   PetscInt  ntsteps; /* Number of time steps */
     32:   PetscInt  nvars;   /* Number of variables in the equation system */
     33:   PetscBool io;
     34: };
     35: typedef pUserCtx *UserCtx;

     37: PetscErrorCode Initialize_AppContext(UserCtx *puser)
     38: {
     39:   UserCtx user;

     41:   PetscFunctionBegin;
     42:   PetscCall(PetscNew(&user));
     43:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "ex35.cxx");
     44:   {
     45:     user->nvars     = 2;
     46:     user->A         = 1;
     47:     user->B         = 3;
     48:     user->alpha     = 0.02;
     49:     user->leftbc.u  = 1;
     50:     user->rightbc.u = 1;
     51:     user->leftbc.v  = 3;
     52:     user->rightbc.v = 3;
     53:     user->n         = 10;
     54:     user->ntsteps   = 10000;
     55:     user->io        = PETSC_FALSE;
     56:     PetscCall(PetscOptionsReal("-A", "Reaction rate", "ex35.cxx", user->A, &user->A, NULL));
     57:     PetscCall(PetscOptionsReal("-B", "Reaction rate", "ex35.cxx", user->B, &user->B, NULL));
     58:     PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "ex35.cxx", user->alpha, &user->alpha, NULL));
     59:     PetscCall(PetscOptionsScalar("-uleft", "Dirichlet boundary condition", "ex35.cxx", user->leftbc.u, &user->leftbc.u, NULL));
     60:     PetscCall(PetscOptionsScalar("-uright", "Dirichlet boundary condition", "ex35.cxx", user->rightbc.u, &user->rightbc.u, NULL));
     61:     PetscCall(PetscOptionsScalar("-vleft", "Dirichlet boundary condition", "ex35.cxx", user->leftbc.v, &user->leftbc.v, NULL));
     62:     PetscCall(PetscOptionsScalar("-vright", "Dirichlet boundary condition", "ex35.cxx", user->rightbc.v, &user->rightbc.v, NULL));
     63:     PetscCall(PetscOptionsInt("-n", "Number of 1-D elements", "ex35.cxx", user->n, &user->n, NULL));
     64:     PetscCall(PetscOptionsInt("-ndt", "Number of time steps", "ex35.cxx", user->ntsteps, &user->ntsteps, NULL));
     65:     PetscCall(PetscOptionsBool("-io", "Write the mesh and solution output to a file.", "ex35.cxx", user->io, &user->io, NULL));
     66:     user->npts = user->n + 1;
     67:   }
     68:   PetscOptionsEnd();

     70:   *puser = user;
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: PetscErrorCode Destroy_AppContext(UserCtx *user)
     75: {
     76:   PetscFunctionBegin;
     77:   PetscCall(PetscFree(*user));
     78:   PetscFunctionReturn(PETSC_SUCCESS);
     79: }

     81: static PetscErrorCode FormInitialSolution(TS, Vec, void *);
     82: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
     83: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
     84: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);

     86: /****************
     87:  *              *
     88:  *     MAIN     *
     89:  *              *
     90:  ****************/
     91: int main(int argc, char **argv)
     92: {
     93:   TS                ts; /* nonlinear solver */
     94:   Vec               X;  /* solution, residual vectors */
     95:   Mat               J;  /* Jacobian matrix */
     96:   PetscInt          steps, mx;
     97:   PetscReal         hx, dt, ftime;
     98:   UserCtx           user; /* user-defined work context */
     99:   TSConvergedReason reason;
    100:   DM                dm;
    101:   const char       *fields[2] = {"U", "V"};

    103:   PetscFunctionBeginUser;
    104:   PetscCall(PetscInitialize(&argc, &argv, nullptr, help));

    106:   /* Initialize the user context struct */
    107:   PetscCall(Initialize_AppContext(&user));

    109:   /* Fill in the user defined work context: */
    110:   PetscCall(DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, NULL, user->n, 1, &dm));
    111:   PetscCall(DMMoabSetFieldNames(dm, user->nvars, fields));
    112:   PetscCall(DMMoabSetBlockSize(dm, user->nvars));
    113:   PetscCall(DMSetFromOptions(dm));

    115:   /* SetUp the data structures for DMMOAB */
    116:   PetscCall(DMSetUp(dm));

    118:   /*  Create timestepping solver context */
    119:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
    120:   PetscCall(TSSetDM(ts, dm));
    121:   PetscCall(TSSetType(ts, TSARKIMEX));
    122:   PetscCall(TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1));
    123:   PetscCall(DMSetMatType(dm, MATBAIJ));
    124:   PetscCall(DMCreateMatrix(dm, &J));

    126:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, user));
    127:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, user));
    128:   PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, user));

    130:   ftime = 10.0;
    131:   PetscCall(TSSetMaxSteps(ts, user->ntsteps));
    132:   PetscCall(TSSetMaxTime(ts, ftime));
    133:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

    135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    136:      Create the solution vector and set the initial conditions
    137:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    138:   PetscCall(DMCreateGlobalVector(dm, &X));

    140:   PetscCall(FormInitialSolution(ts, X, user));
    141:   PetscCall(TSSetSolution(ts, X));
    142:   PetscCall(VecGetSize(X, &mx));
    143:   hx = 1.0 / (PetscReal)(mx / 2 - 1);
    144:   dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */
    145:   PetscCall(TSSetTimeStep(ts, dt));

    147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    148:      Set runtime options
    149:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    150:   PetscCall(TSSetFromOptions(ts));

    152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    153:      Solve nonlinear system
    154:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    155:   PetscCall(TSSolve(ts, X));
    156:   PetscCall(TSGetSolveTime(ts, &ftime));
    157:   PetscCall(TSGetStepNumber(ts, &steps));
    158:   PetscCall(TSGetConvergedReason(ts, &reason));
    159:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

    161:   if (user->io) {
    162:     /* Print the numerical solution to screen and then dump to file */
    163:     PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));

    165:     /* Write out the solution along with the mesh */
    166:     PetscCall(DMMoabSetGlobalFieldVector(dm, X));
    167: #ifdef MOAB_HAVE_HDF5
    168:     PetscCall(DMMoabOutput(dm, "ex35.h5m", ""));
    169: #else
    170:     /* MOAB does not support true parallel writers that aren't HDF5 based
    171:        And so if you are using VTK as the output format in parallel,
    172:        the data could be jumbled due to the order in which the processors
    173:        write out their parts of the mesh and solution tags
    174:     */
    175:     PetscCall(DMMoabOutput(dm, "ex35.vtk", ""));
    176: #endif
    177:   }

    179:   /* Free work space.
    180:      Free all PETSc related resources: */
    181:   PetscCall(MatDestroy(&J));
    182:   PetscCall(VecDestroy(&X));
    183:   PetscCall(TSDestroy(&ts));
    184:   PetscCall(DMDestroy(&dm));

    186:   /* Free all MOAB related resources: */
    187:   PetscCall(Destroy_AppContext(&user));
    188:   PetscCall(PetscFinalize());
    189:   return 0;
    190: }

    192: /*
    193:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
    194: */
    195: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
    196: {
    197:   UserCtx            user = (UserCtx)ptr;
    198:   PetscInt           dof;
    199:   PetscReal          hx;
    200:   DM                 dm;
    201:   const moab::Range *vlocal;
    202:   PetscBool          vonboundary;

    204:   PetscFunctionBegin;
    205:   PetscCall(TSGetDM(ts, &dm));

    207:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
    208:   PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL));

    210:   /* compute local element sizes - structured grid */
    211:   hx = 1.0 / user->n;

    213:   /* Compute function over the locally owned part of the grid
    214:      Assemble the operator by looping over edges and computing
    215:      contribution for each vertex dof                         */
    216:   for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
    217:     const moab::EntityHandle vhandle = *iter;

    219:     PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof));

    221:     /* check if vertex is on the boundary */
    222:     PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &vonboundary));

    224:     if (vonboundary) {
    225:       const PetscScalar bcvals[2][2] = {
    226:         {hx, 0 },
    227:         {0,  hx}
    228:       };
    229:       PetscCall(MatSetValuesBlocked(Jpre, 1, &dof, 1, &dof, &bcvals[0][0], INSERT_VALUES));
    230:     } else {
    231:       const PetscInt    row = dof, col[] = {dof - 1, dof, dof + 1};
    232:       const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
    233:       const PetscScalar vals[2][3][2] = {
    234:         {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
    235:         {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
    236:       };
    237:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES));
    238:     }
    239:   }

    241:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    242:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    243:   if (J != Jpre) {
    244:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    245:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    246:   }
    247:   PetscFunctionReturn(PETSC_SUCCESS);
    248: }

    250: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
    251: {
    252:   UserCtx            user = (UserCtx)ptr;
    253:   DM                 dm;
    254:   PetscReal          hx;
    255:   const Field       *x;
    256:   Field             *f;
    257:   PetscInt           dof;
    258:   const moab::Range *ownedvtx;

    260:   PetscFunctionBegin;
    261:   hx = 1.0 / user->n;
    262:   PetscCall(TSGetDM(ts, &dm));

    264:   /* Get pointers to vector data */
    265:   PetscCall(VecSet(F, 0.0));

    267:   PetscCall(DMMoabVecGetArrayRead(dm, X, &x));
    268:   PetscCall(DMMoabVecGetArray(dm, F, &f));

    270:   PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));

    272:   /* Compute function over the locally owned part of the grid */
    273:   for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
    274:     const moab::EntityHandle vhandle = *iter;
    275:     PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof));

    277:     PetscScalar u = x[dof].u, v = x[dof].v;
    278:     f[dof].u = hx * (user->A + u * u * v - (user->B + 1) * u);
    279:     f[dof].v = hx * (user->B * u - u * u * v);
    280:   }

    282:   /* Restore vectors */
    283:   PetscCall(DMMoabVecRestoreArrayRead(dm, X, &x));
    284:   PetscCall(DMMoabVecRestoreArray(dm, F, &f));
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
    289: {
    290:   UserCtx            user = (UserCtx)ctx;
    291:   DM                 dm;
    292:   Field             *x, *xdot, *f;
    293:   PetscReal          hx;
    294:   Vec                Xloc;
    295:   PetscInt           i, bcindx;
    296:   PetscBool          elem_on_boundary;
    297:   const moab::Range *vlocal;

    299:   PetscFunctionBegin;
    300:   hx = 1.0 / user->n;
    301:   PetscCall(TSGetDM(ts, &dm));

    303:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
    304:   PetscCall(DMMoabGetLocalVertices(dm, &vlocal, NULL));

    306:   /* reset the residual vector */
    307:   PetscCall(VecSet(F, 0.0));

    309:   PetscCall(DMGetLocalVector(dm, &Xloc));
    310:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
    311:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));

    313:   /* get the local representation of the arrays from Vectors */
    314:   PetscCall(DMMoabVecGetArrayRead(dm, Xloc, &x));
    315:   PetscCall(DMMoabVecGetArrayRead(dm, Xdot, &xdot));
    316:   PetscCall(DMMoabVecGetArray(dm, F, &f));

    318:   /* loop over local elements */
    319:   for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
    320:     const moab::EntityHandle vhandle = *iter;

    322:     PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &i));

    324:     /* check if vertex is on the boundary */
    325:     PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &elem_on_boundary));

    327:     if (elem_on_boundary) {
    328:       PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx));
    329:       if (bcindx == 0) { /* Apply left BC */
    330:         f[i].u = hx * (x[i].u - user->leftbc.u);
    331:         f[i].v = hx * (x[i].v - user->leftbc.v);
    332:       } else { /* Apply right BC */
    333:         f[i].u = hx * (x[i].u - user->rightbc.u);
    334:         f[i].v = hx * (x[i].v - user->rightbc.v);
    335:       }
    336:     } else {
    337:       f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
    338:       f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
    339:     }
    340:   }

    342:   /* Restore data */
    343:   PetscCall(DMMoabVecRestoreArrayRead(dm, Xloc, &x));
    344:   PetscCall(DMMoabVecRestoreArrayRead(dm, Xdot, &xdot));
    345:   PetscCall(DMMoabVecRestoreArray(dm, F, &f));
    346:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
    347:   PetscFunctionReturn(PETSC_SUCCESS);
    348: }

    350: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
    351: {
    352:   UserCtx               user = (UserCtx)ctx;
    353:   PetscReal             vpos[3];
    354:   DM                    dm;
    355:   Field                *x;
    356:   const moab::Range    *vowned;
    357:   PetscInt              dof;
    358:   moab::Range::iterator iter;

    360:   PetscFunctionBegin;
    361:   PetscCall(TSGetDM(ts, &dm));

    363:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
    364:   PetscCall(DMMoabGetLocalVertices(dm, &vowned, NULL));

    366:   PetscCall(VecSet(X, 0.0));

    368:   /* Get pointers to vector data */
    369:   PetscCall(DMMoabVecGetArray(dm, X, &x));

    371:   /* Compute function over the locally owned part of the grid */
    372:   for (moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) {
    373:     const moab::EntityHandle vhandle = *iter;
    374:     PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof));

    376:     /* compute the mid-point of the element and use a 1-point lumped quadrature */
    377:     PetscCall(DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos));

    379:     PetscReal xi = vpos[0];
    380:     x[dof].u     = user->leftbc.u * (1. - xi) + user->rightbc.u * xi + PetscSinReal(2. * PETSC_PI * xi);
    381:     x[dof].v     = user->leftbc.v * (1. - xi) + user->rightbc.v * xi;
    382:   }

    384:   /* Restore vectors */
    385:   PetscCall(DMMoabVecRestoreArray(dm, X, &x));
    386:   PetscFunctionReturn(PETSC_SUCCESS);
    387: }

    389: /*TEST

    391:     build:
    392:       requires: moab !complex

    394:     test:
    395:       args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none

    397:     test:
    398:       suffix: 2
    399:       nsize: 2
    400:       args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io
    401:       TODO:

    403: TEST*/