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