Actual source code: dmlocalts.c

      1: #include <petsc/private/dmimpl.h>
      2: #include <petsc/private/tsimpl.h>

      4: typedef struct {
      5:   PetscErrorCode (*boundarylocal)(DM, PetscReal, Vec, Vec, void *);
      6:   PetscErrorCode (*ifunctionlocal)(DM, PetscReal, Vec, Vec, Vec, void *);
      7:   PetscErrorCode (*ijacobianlocal)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
      8:   PetscErrorCode (*rhsfunctionlocal)(DM, PetscReal, Vec, Vec, void *);
      9:   void *boundarylocalctx;
     10:   void *ifunctionlocalctx;
     11:   void *ijacobianlocalctx;
     12:   void *rhsfunctionlocalctx;
     13:   Vec   lumpedmassinv;
     14:   Mat   mass;
     15:   KSP   kspmass;
     16: } DMTS_Local;

     18: static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm)
     19: {
     20:   PetscFunctionBegin;
     21:   PetscCall(PetscFree(tdm->data));
     22:   PetscFunctionReturn(PETSC_SUCCESS);
     23: }

     25: static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
     26: {
     27:   PetscFunctionBegin;
     28:   PetscCall(PetscNew((DMTS_Local **)&tdm->data));
     29:   if (oldtdm->data) PetscCall(PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local)));
     30:   PetscFunctionReturn(PETSC_SUCCESS);
     31: }

     33: static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts)
     34: {
     35:   PetscFunctionBegin;
     36:   *dmlocalts = NULL;
     37:   if (!tdm->data) {
     38:     PetscCall(PetscNew((DMTS_Local **)&tdm->data));

     40:     tdm->ops->destroy   = DMTSDestroy_DMLocal;
     41:     tdm->ops->duplicate = DMTSDuplicate_DMLocal;
     42:   }
     43:   *dmlocalts = (DMTS_Local *)tdm->data;
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, void *ctx)
     48: {
     49:   DM          dm;
     50:   Vec         locX, locX_t, locF;
     51:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

     53:   PetscFunctionBegin;
     58:   PetscCall(TSGetDM(ts, &dm));
     59:   PetscCall(DMGetLocalVector(dm, &locX));
     60:   PetscCall(DMGetLocalVector(dm, &locX_t));
     61:   PetscCall(DMGetLocalVector(dm, &locF));
     62:   PetscCall(VecZeroEntries(locX));
     63:   PetscCall(VecZeroEntries(locX_t));
     64:   if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx));
     65:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
     66:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
     67:   PetscCall(DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t));
     68:   PetscCall(DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t));
     69:   PetscCall(VecZeroEntries(locF));
     70:   CHKMEMQ;
     71:   PetscCall((*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx));
     72:   CHKMEMQ;
     73:   PetscCall(VecZeroEntries(F));
     74:   PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F));
     75:   PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F));
     76:   PetscCall(DMRestoreLocalVector(dm, &locX));
     77:   PetscCall(DMRestoreLocalVector(dm, &locX_t));
     78:   PetscCall(DMRestoreLocalVector(dm, &locF));
     79:   PetscFunctionReturn(PETSC_SUCCESS);
     80: }

     82: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, void *ctx)
     83: {
     84:   DM          dm;
     85:   Vec         locX, locF;
     86:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

     88:   PetscFunctionBegin;
     92:   PetscCall(TSGetDM(ts, &dm));
     93:   PetscCall(DMGetLocalVector(dm, &locX));
     94:   PetscCall(DMGetLocalVector(dm, &locF));
     95:   PetscCall(VecZeroEntries(locX));
     96:   if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, NULL, dmlocalts->boundarylocalctx));
     97:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
     98:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
     99:   PetscCall(VecZeroEntries(locF));
    100:   CHKMEMQ;
    101:   PetscCall((*dmlocalts->rhsfunctionlocal)(dm, time, locX, locF, dmlocalts->rhsfunctionlocalctx));
    102:   CHKMEMQ;
    103:   PetscCall(VecZeroEntries(F));
    104:   PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F));
    105:   PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F));
    106:   if (dmlocalts->lumpedmassinv) {
    107:     PetscCall(VecPointwiseMult(F, dmlocalts->lumpedmassinv, F));
    108:   } else if (dmlocalts->kspmass) {
    109:     Vec tmp;

    111:     PetscCall(DMGetGlobalVector(dm, &tmp));
    112:     PetscCall(KSPSolve(dmlocalts->kspmass, F, tmp));
    113:     PetscCall(VecCopy(tmp, F));
    114:     PetscCall(DMRestoreGlobalVector(dm, &tmp));
    115:   }
    116:   PetscCall(DMRestoreLocalVector(dm, &locX));
    117:   PetscCall(DMRestoreLocalVector(dm, &locF));
    118:   PetscFunctionReturn(PETSC_SUCCESS);
    119: }

    121: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, void *ctx)
    122: {
    123:   DM          dm;
    124:   Vec         locX, locX_t;
    125:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

    127:   PetscFunctionBegin;
    128:   PetscCall(TSGetDM(ts, &dm));
    129:   if (dmlocalts->ijacobianlocal) {
    130:     PetscCall(DMGetLocalVector(dm, &locX));
    131:     PetscCall(DMGetLocalVector(dm, &locX_t));
    132:     PetscCall(VecZeroEntries(locX));
    133:     PetscCall(VecZeroEntries(locX_t));
    134:     if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx));
    135:     PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
    136:     PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
    137:     PetscCall(DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t));
    138:     PetscCall(DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t));
    139:     CHKMEMQ;
    140:     PetscCall((*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx));
    141:     CHKMEMQ;
    142:     PetscCall(DMRestoreLocalVector(dm, &locX));
    143:     PetscCall(DMRestoreLocalVector(dm, &locX_t));
    144:   } else {
    145:     MatFDColoring fdcoloring;
    146:     PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring));
    147:     if (!fdcoloring) {
    148:       ISColoring coloring;

    150:       PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
    151:       PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
    152:       PetscCall(ISColoringDestroy(&coloring));
    153:       switch (dm->coloringtype) {
    154:       case IS_COLORING_GLOBAL:
    155:         PetscCall(MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))TSComputeIFunction_DMLocal, dmlocalts));
    156:         break;
    157:       default:
    158:         SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
    159:       }
    160:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
    161:       PetscCall(MatFDColoringSetFromOptions(fdcoloring));
    162:       PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
    163:       PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
    164:       PetscCall(PetscObjectDereference((PetscObject)fdcoloring));

    166:       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
    167:        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
    168:        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
    169:        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
    170:        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
    171:        */
    172:       PetscCall(PetscObjectDereference((PetscObject)dm));
    173:     }
    174:     PetscCall(MatFDColoringApply(B, fdcoloring, X, ts));
    175:   }
    176:   /* This will be redundant if the user called both, but it's too common to forget. */
    177:   if (A != B) {
    178:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
    179:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
    180:   }
    181:   PetscFunctionReturn(PETSC_SUCCESS);
    182: }

    184: /*@C
    185:   DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.

    187:   Logically Collective

    189:   Input Parameters:
    190: + dm   - `DM` to associate callback with
    191: . func - local function evaluation
    192: - ctx  - context for function evaluation

    194:   Level: intermediate

    196:   Notes:
    197:   `func` should set the essential boundary data for the local portion of the solution, as
    198:   well its time derivative (if it is not `NULL`).

    200:   Vectors are initialized to zero before this function, so it is only needed for non
    201:   homogeneous data.

    203:   This function is somewhat optional: boundary data could potentially be inserted by a function
    204:   passed to `DMTSSetIFunctionLocal()`. The use case for this function is for discretizations
    205:   with constraints (see `DMGetDefaultConstraints()`): this function inserts boundary values
    206:   before constraint interpolation.

    208: .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    209: @*/
    210: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
    211: {
    212:   DMTS        tdm;
    213:   DMTS_Local *dmlocalts;

    215:   PetscFunctionBegin;
    217:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    218:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    220:   dmlocalts->boundarylocal    = func;
    221:   dmlocalts->boundarylocalctx = ctx;
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

    225: /*@C
    226:   DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
    227:   containing the local vector information PLUS ghost point information. It should compute a result for all local
    228:   elements and `DM` will automatically accumulate the overlapping values.

    230:   Logically Collective

    232:   Input Parameter:
    233: . dm - `DM` to associate callback with

    235:   Output Parameters:
    236: + func - local function evaluation
    237: - ctx  - context for function evaluation

    239:   Level: beginner

    241: .seealso: [](ch_ts), `DM`, `DMTSSetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    242: @*/
    243: PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx)
    244: {
    245:   DMTS        tdm;
    246:   DMTS_Local *dmlocalts;

    248:   PetscFunctionBegin;
    250:   PetscCall(DMGetDMTS(dm, &tdm));
    251:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    252:   if (func) {
    253:     PetscAssertPointer(func, 2);
    254:     *func = dmlocalts->ifunctionlocal;
    255:   }
    256:   if (ctx) {
    257:     PetscAssertPointer(ctx, 3);
    258:     *ctx = dmlocalts->ifunctionlocalctx;
    259:   }
    260:   PetscFunctionReturn(PETSC_SUCCESS);
    261: }

    263: /*@C
    264:   DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
    265:   containing the local vector information PLUS ghost point information. It should compute a result for all local
    266:   elements and `DM` will automatically accumulate the overlapping values.

    268:   Logically Collective

    270:   Input Parameters:
    271: + dm   - `DM` to associate callback with
    272: . func - local function evaluation
    273: - ctx  - context for function evaluation

    275:   Level: beginner

    277: .seealso: [](ch_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    278: @*/
    279: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
    280: {
    281:   DMTS        tdm;
    282:   DMTS_Local *dmlocalts;

    284:   PetscFunctionBegin;
    286:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    287:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    289:   dmlocalts->ifunctionlocal    = func;
    290:   dmlocalts->ifunctionlocalctx = ctx;

    292:   PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts));
    293:   if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
    294:     PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
    295:   }
    296:   PetscFunctionReturn(PETSC_SUCCESS);
    297: }

    299: /*@C
    300:   DMTSGetIJacobianLocal - get a local Jacobian evaluation function

    302:   Logically Collective

    304:   Input Parameter:
    305: . dm - `DM` to associate callback with

    307:   Output Parameters:
    308: + func - local Jacobian evaluation
    309: - ctx  - optional context for local Jacobian evaluation

    311:   Level: beginner

    313: .seealso: [](ch_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
    314: @*/
    315: PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void **ctx)
    316: {
    317:   DMTS        tdm;
    318:   DMTS_Local *dmlocalts;

    320:   PetscFunctionBegin;
    322:   PetscCall(DMGetDMTS(dm, &tdm));
    323:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    324:   if (func) {
    325:     PetscAssertPointer(func, 2);
    326:     *func = dmlocalts->ijacobianlocal;
    327:   }
    328:   if (ctx) {
    329:     PetscAssertPointer(ctx, 3);
    330:     *ctx = dmlocalts->ijacobianlocalctx;
    331:   }
    332:   PetscFunctionReturn(PETSC_SUCCESS);
    333: }

    335: /*@C
    336:   DMTSSetIJacobianLocal - set a local Jacobian evaluation function

    338:   Logically Collective

    340:   Input Parameters:
    341: + dm   - `DM` to associate callback with
    342: . func - local Jacobian evaluation
    343: - ctx  - optional context for local Jacobian evaluation

    345:   Level: beginner

    347: .seealso: [](ch_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
    348: @*/
    349: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
    350: {
    351:   DMTS        tdm;
    352:   DMTS_Local *dmlocalts;

    354:   PetscFunctionBegin;
    356:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    357:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    359:   dmlocalts->ijacobianlocal    = func;
    360:   dmlocalts->ijacobianlocalctx = ctx;

    362:   PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
    363:   PetscFunctionReturn(PETSC_SUCCESS);
    364: }

    366: /*@C
    367:   DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
    368:   containing the local vector information PLUS ghost point information. It should compute a result for all local
    369:   elements and `DM` will automatically accumulate the overlapping values.

    371:   Logically Collective

    373:   Input Parameter:
    374: . dm - `DM` to associate callback with

    376:   Output Parameters:
    377: + func - local function evaluation
    378: - ctx  - context for function evaluation

    380:   Level: beginner

    382: .seealso: [](ch_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    383: @*/
    384: PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx)
    385: {
    386:   DMTS        tdm;
    387:   DMTS_Local *dmlocalts;

    389:   PetscFunctionBegin;
    391:   PetscCall(DMGetDMTS(dm, &tdm));
    392:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    393:   if (func) {
    394:     PetscAssertPointer(func, 2);
    395:     *func = dmlocalts->rhsfunctionlocal;
    396:   }
    397:   if (ctx) {
    398:     PetscAssertPointer(ctx, 3);
    399:     *ctx = dmlocalts->rhsfunctionlocalctx;
    400:   }
    401:   PetscFunctionReturn(PETSC_SUCCESS);
    402: }

    404: /*@C
    405:   DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
    406:   containing the local vector information PLUS ghost point information. It should compute a result for all local
    407:   elements and `DM` will automatically accumulate the overlapping values.

    409:   Logically Collective

    411:   Input Parameters:
    412: + dm   - `DM` to associate callback with
    413: . func - local function evaluation
    414: - ctx  - context for function evaluation

    416:   Level: beginner

    418: .seealso: [](ch_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    419: @*/
    420: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
    421: {
    422:   DMTS        tdm;
    423:   DMTS_Local *dmlocalts;

    425:   PetscFunctionBegin;
    427:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    428:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    430:   dmlocalts->rhsfunctionlocal    = func;
    431:   dmlocalts->rhsfunctionlocalctx = ctx;

    433:   PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts));
    434:   PetscFunctionReturn(PETSC_SUCCESS);
    435: }

    437: /*@
    438:   DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.

    440:   Collective

    442:   Input Parameter:
    443: . dm - `DM` providing the mass matrix

    445:   Level: developer

    447:   Note:
    448:   The idea here is that an explicit system can be given a mass matrix, based on the `DM`, which is inverted on the RHS at each step.

    450: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
    451: @*/
    452: PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
    453: {
    454:   DMTS        tdm;
    455:   DMTS_Local *dmlocalts;
    456:   const char *prefix;

    458:   PetscFunctionBegin;
    460:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    461:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    462:   PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass));
    463:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass));
    464:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    465:   PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix));
    466:   PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_"));
    467:   PetscCall(KSPSetFromOptions(dmlocalts->kspmass));
    468:   PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass));
    469:   PetscFunctionReturn(PETSC_SUCCESS);
    470: }

    472: /*@
    473:   DMTSCreateRHSMassMatrixLumped - This creates the lumped mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.

    475:   Collective

    477:   Input Parameter:
    478: . dm - `DM` providing the mass matrix

    480:   Level: developer

    482:   Note:
    483:   The idea here is that an explicit system can be given a mass matrix, based on the `DM`, which is inverted on the RHS at each step.
    484:   Since the matrix is lumped, inversion is trivial.

    486: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
    487: @*/
    488: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
    489: {
    490:   DMTS        tdm;
    491:   DMTS_Local *dmlocalts;

    493:   PetscFunctionBegin;
    495:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    496:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    497:   PetscCall(DMCreateMassMatrixLumped(dm, NULL, &dmlocalts->lumpedmassinv));
    498:   PetscCall(VecReciprocal(dmlocalts->lumpedmassinv));
    499:   PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view"));
    500:   PetscFunctionReturn(PETSC_SUCCESS);
    501: }

    503: /*@
    504:   DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist.

    506:   Logically Collective

    508:   Input Parameter:
    509: . dm - `DM` providing the mass matrix

    511:   Level: developer

    513: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMTS`
    514: @*/
    515: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
    516: {
    517:   DMTS        tdm;
    518:   DMTS_Local *dmlocalts;

    520:   PetscFunctionBegin;
    522:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    523:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    524:   PetscCall(VecDestroy(&dmlocalts->lumpedmassinv));
    525:   PetscCall(MatDestroy(&dmlocalts->mass));
    526:   PetscCall(KSPDestroy(&dmlocalts->kspmass));
    527:   PetscFunctionReturn(PETSC_SUCCESS);
    528: }