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

     80:   /* remove nullspace from residual */
     81:   {
     82:     MatNullSpace nullsp;
     83:     PetscCall(PetscObjectQuery((PetscObject)dm, "__dmtsnullspace", (PetscObject *)&nullsp));
     84:     if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, F));
     85:   }
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, PetscCtx ctx)
     90: {
     91:   DM          dm;
     92:   Vec         locX, locF;
     93:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

     95:   PetscFunctionBegin;
     99:   PetscCall(TSGetDM(ts, &dm));
    100:   PetscCall(DMGetLocalVector(dm, &locX));
    101:   PetscCall(DMGetLocalVector(dm, &locF));
    102:   PetscCall(VecZeroEntries(locX));
    103:   if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, NULL, dmlocalts->boundarylocalctx));
    104:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
    105:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
    106:   PetscCall(VecZeroEntries(locF));
    107:   CHKMEMQ;
    108:   PetscCall((*dmlocalts->rhsfunctionlocal)(dm, time, locX, locF, dmlocalts->rhsfunctionlocalctx));
    109:   CHKMEMQ;
    110:   PetscCall(VecZeroEntries(F));
    111:   PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F));
    112:   PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F));
    113:   if (dmlocalts->lumpedmassinv) PetscCall(VecPointwiseMult(F, dmlocalts->lumpedmassinv, F));
    114:   else if (dmlocalts->kspmass) {
    115:     Vec tmp;

    117:     PetscCall(DMGetGlobalVector(dm, &tmp));
    118:     PetscCall(KSPSolve(dmlocalts->kspmass, F, tmp));
    119:     PetscCall(VecCopy(tmp, F));
    120:     PetscCall(DMRestoreGlobalVector(dm, &tmp));
    121:   }
    122:   PetscCall(DMRestoreLocalVector(dm, &locX));
    123:   PetscCall(DMRestoreLocalVector(dm, &locF));
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, PetscCtx ctx)
    128: {
    129:   DM          dm;
    130:   Vec         locX, locX_t;
    131:   DMTS_Local *dmlocalts = (DMTS_Local *)ctx;

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

    156:       PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
    157:       PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
    158:       PetscCall(ISColoringDestroy(&coloring));
    159:       switch (dm->coloringtype) {
    160:       case IS_COLORING_GLOBAL:
    161:         PetscCall(MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn *)(PetscVoidFn *)TSComputeIFunction_DMLocal, dmlocalts));
    162:         break;
    163:       default:
    164:         SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
    165:       }
    166:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
    167:       PetscCall(MatFDColoringSetFromOptions(fdcoloring));
    168:       PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
    169:       PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
    170:       PetscCall(PetscObjectDereference((PetscObject)fdcoloring));

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

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

    193:   Logically Collective

    195:   Input Parameters:
    196: + dm   - `DM` to associate callback with
    197: . func - local function evaluation
    198: - ctx  - context for function evaluation

    200:   Calling sequence of `func`:
    201: + dm  - the `DM`
    202: . t   - the current time
    203: . u   - the current solution
    204: . f   - output, the computed right hand side function
    205: - ctx - the application context for the function

    207:   Level: intermediate

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

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

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

    221: .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    222: @*/
    223: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec f, PetscCtx ctx), PetscCtx ctx)
    224: {
    225:   DMTS        tdm;
    226:   DMTS_Local *dmlocalts;

    228:   PetscFunctionBegin;
    230:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    231:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    233:   dmlocalts->boundarylocal    = func;
    234:   dmlocalts->boundarylocalctx = ctx;
    235:   PetscFunctionReturn(PETSC_SUCCESS);
    236: }

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

    243:   Logically Collective

    245:   Input Parameter:
    246: . dm - `DM` to associate callback with

    248:   Output Parameters:
    249: + func - local function evaluation
    250: - ctx  - context for function evaluation

    252:   Calling sequence of `func`:
    253: + dm   - the `DM`
    254: . t    - the current time
    255: . u    - the current solution
    256: . udot - the derivative of `u`
    257: . F    - output, the computed implicit function
    258: - ctx  - the application context for the function

    260:   Level: beginner

    262: .seealso: [](ch_ts), `DM`, `DMTSSetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    263: @*/
    264: PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM dm, PetscReal t, Vec u, Vec udot, Vec F, PetscCtx ctx), PetscCtxRt ctx)
    265: {
    266:   DMTS        tdm;
    267:   DMTS_Local *dmlocalts;

    269:   PetscFunctionBegin;
    271:   PetscCall(DMGetDMTS(dm, &tdm));
    272:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    273:   if (func) {
    274:     PetscAssertPointer(func, 2);
    275:     *func = dmlocalts->ifunctionlocal;
    276:   }
    277:   if (ctx) {
    278:     PetscAssertPointer(ctx, 3);
    279:     *(void **)ctx = dmlocalts->ifunctionlocalctx;
    280:   }
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }

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

    289:   Logically Collective

    291:   Input Parameters:
    292: + dm   - `DM` to associate callback with
    293: . func - local function evaluation
    294: - ctx  - context for function evaluation

    296:   Calling sequence of `func`:
    297: + dm   - the `DM`
    298: . t    - the current time
    299: . u    - the current solution
    300: . udot - the derivative of `u`
    301: . F    - output, the computed implicit function
    302: - ctx  - the application context for the function

    304:   Level: beginner

    306: .seealso: [](ch_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    307: @*/
    308: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec udot, Vec F, PetscCtx ctx), PetscCtx ctx)
    309: {
    310:   DMTS        tdm;
    311:   DMTS_Local *dmlocalts;

    313:   PetscFunctionBegin;
    315:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    316:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    318:   dmlocalts->ifunctionlocal    = func;
    319:   dmlocalts->ifunctionlocalctx = ctx;

    321:   PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts));
    322:   if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
    323:     PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
    324:   }
    325:   PetscFunctionReturn(PETSC_SUCCESS);
    326: }

    328: /*@C
    329:   DMTSGetIJacobianLocal - get a local Jacobian evaluation function

    331:   Logically Collective

    333:   Input Parameter:
    334: . dm - `DM` to associate callback with

    336:   Output Parameters:
    337: + func - local Jacobian evaluation
    338: - ctx  - optional context for local Jacobian evaluation

    340:   Calling sequence of `func`:
    341: + dm    - the `DM`
    342: . t     - the current time
    343: . u     - the current solution
    344: . udot  - the derivative of `u`
    345: . shift - the shift factoring arising from the implicit time-step
    346: . J     - output, the Jacobian
    347: . Jpre  - output, matrix from which to compute the preconditioner for `J`, often the same as `J`
    348: - ctx   - the application context for the function

    350:   Level: beginner

    352: .seealso: [](ch_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
    353: @*/
    354: PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM dm, PetscReal t, Vec u, Vec udot, PetscReal shift, Mat J, Mat Jpre, PetscCtx ctx), PetscCtxRt ctx)
    355: {
    356:   DMTS        tdm;
    357:   DMTS_Local *dmlocalts;

    359:   PetscFunctionBegin;
    361:   PetscCall(DMGetDMTS(dm, &tdm));
    362:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    363:   if (func) {
    364:     PetscAssertPointer(func, 2);
    365:     *func = dmlocalts->ijacobianlocal;
    366:   }
    367:   if (ctx) {
    368:     PetscAssertPointer(ctx, 3);
    369:     *(void **)ctx = dmlocalts->ijacobianlocalctx;
    370:   }
    371:   PetscFunctionReturn(PETSC_SUCCESS);
    372: }

    374: /*@C
    375:   DMTSSetIJacobianLocal - set a local Jacobian evaluation function

    377:   Logically Collective

    379:   Input Parameters:
    380: + dm   - `DM` to associate callback with
    381: . func - local Jacobian evaluation
    382: - ctx  - optional context for local Jacobian evaluation

    384:   Calling sequence of `func`:
    385: + dm    - the `DM`
    386: . t     - the current time
    387: . u     - the current solution
    388: . udot  - the derivative of `u`
    389: . shift - the shift factoring arising from the implicit time-step
    390: . J     - output, the Jacobian
    391: . Jpre  - output, matrix from which to compute the preconditioner for `J`, often the same as `J`
    392: - ctx   - the application context for the function

    394:   Level: beginner

    396: .seealso: [](ch_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
    397: @*/
    398: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec udot, PetscReal shift, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
    399: {
    400:   DMTS        tdm;
    401:   DMTS_Local *dmlocalts;

    403:   PetscFunctionBegin;
    405:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    406:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    408:   dmlocalts->ijacobianlocal    = func;
    409:   dmlocalts->ijacobianlocalctx = ctx;

    411:   PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
    412:   PetscFunctionReturn(PETSC_SUCCESS);
    413: }

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

    420:   Logically Collective

    422:   Input Parameter:
    423: . dm - `DM` to associate callback with

    425:   Output Parameters:
    426: + func - local function evaluation
    427: - ctx  - context for function evaluation

    429:   Calling sequence of `func`:
    430: + dm   - the `DM`
    431: . t    - the current time
    432: . u    - the current solution
    433: . udot - output, the evaluated right hand side
    434: - ctx  - the application context for the function

    436:   Level: beginner

    438: .seealso: [](ch_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    439: @*/
    440: PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM dm, PetscReal t, Vec u, Vec udot, PetscCtx ctx), PetscCtxRt ctx)
    441: {
    442:   DMTS        tdm;
    443:   DMTS_Local *dmlocalts;

    445:   PetscFunctionBegin;
    447:   PetscCall(DMGetDMTS(dm, &tdm));
    448:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    449:   if (func) {
    450:     PetscAssertPointer(func, 2);
    451:     *func = dmlocalts->rhsfunctionlocal;
    452:   }
    453:   if (ctx) {
    454:     PetscAssertPointer(ctx, 3);
    455:     *(void **)ctx = dmlocalts->rhsfunctionlocalctx;
    456:   }
    457:   PetscFunctionReturn(PETSC_SUCCESS);
    458: }

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

    465:   Logically Collective

    467:   Input Parameters:
    468: + dm   - `DM` to associate callback with
    469: . func - local function evaluation
    470: - ctx  - context for function evaluation

    472:   Calling sequence of `func`:
    473: + dm  - the `DM`
    474: . t   - the current time
    475: . u   - the current solution
    476: . f   - output, the evaluated right hand side
    477: - ctx - the application context for the function

    479:   Level: beginner

    481: .seealso: [](ch_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
    482: @*/
    483: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec f, PetscCtx ctx), PetscCtx ctx)
    484: {
    485:   DMTS        tdm;
    486:   DMTS_Local *dmlocalts;

    488:   PetscFunctionBegin;
    490:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    491:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

    493:   dmlocalts->rhsfunctionlocal    = func;
    494:   dmlocalts->rhsfunctionlocalctx = ctx;

    496:   PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts));
    497:   PetscFunctionReturn(PETSC_SUCCESS);
    498: }

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

    503:   Collective

    505:   Input Parameter:
    506: . dm - `DM` providing the mass matrix

    508:   Level: developer

    510:   Note:
    511:   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.

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

    521:   PetscFunctionBegin;
    523:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    524:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    525:   PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass));
    526:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass));
    527:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    528:   PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix));
    529:   PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_"));
    530:   PetscCall(KSPSetFromOptions(dmlocalts->kspmass));
    531:   PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass));
    532:   PetscFunctionReturn(PETSC_SUCCESS);
    533: }

    535: /*@
    536:   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.

    538:   Collective

    540:   Input Parameter:
    541: . dm - `DM` providing the mass matrix

    543:   Level: developer

    545:   Note:
    546:   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.
    547:   Since the matrix is lumped, inversion is trivial.

    549: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
    550: @*/
    551: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
    552: {
    553:   DMTS        tdm;
    554:   DMTS_Local *dmlocalts;

    556:   PetscFunctionBegin;
    558:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    559:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    560:   PetscCall(DMCreateMassMatrixLumped(dm, NULL, &dmlocalts->lumpedmassinv));
    561:   PetscCall(VecReciprocal(dmlocalts->lumpedmassinv));
    562:   PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view"));
    563:   PetscFunctionReturn(PETSC_SUCCESS);
    564: }

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

    569:   Logically Collective

    571:   Input Parameter:
    572: . dm - `DM` providing the mass matrix

    574:   Level: developer

    576: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMTS`
    577: @*/
    578: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
    579: {
    580:   DMTS        tdm;
    581:   DMTS_Local *dmlocalts;

    583:   PetscFunctionBegin;
    585:   PetscCall(DMGetDMTSWrite(dm, &tdm));
    586:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
    587:   PetscCall(VecDestroy(&dmlocalts->lumpedmassinv));
    588:   PetscCall(MatDestroy(&dmlocalts->mass));
    589:   PetscCall(KSPDestroy(&dmlocalts->kspmass));
    590:   PetscFunctionReturn(PETSC_SUCCESS);
    591: }