xref: /petsc/src/ts/utils/dmlocalts.c (revision fbf9dbe564678ed6eff1806adbc4c4f01b9743f4)
1 #include <petsc/private/dmimpl.h>
2 #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
3 
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;
17 
18 static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm)
19 {
20   PetscFunctionBegin;
21   PetscCall(PetscFree(tdm->data));
22   PetscFunctionReturn(PETSC_SUCCESS);
23 }
24 
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 }
32 
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));
39 
40     tdm->ops->destroy   = DMTSDestroy_DMLocal;
41     tdm->ops->duplicate = DMTSDuplicate_DMLocal;
42   }
43   *dmlocalts = (DMTS_Local *)tdm->data;
44   PetscFunctionReturn(PETSC_SUCCESS);
45 }
46 
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;
52 
53   PetscFunctionBegin;
54   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
55   PetscValidHeaderSpecific(X, VEC_CLASSID, 3);
56   PetscValidHeaderSpecific(X_t, VEC_CLASSID, 4);
57   PetscValidHeaderSpecific(F, VEC_CLASSID, 5);
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 }
81 
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;
87 
88   PetscFunctionBegin;
89   PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
90   PetscValidHeaderSpecific(X, VEC_CLASSID, 3);
91   PetscValidHeaderSpecific(F, VEC_CLASSID, 4);
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;
110 
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 }
120 
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;
126 
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;
149 
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));
165 
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 }
183 
184 /*@C
185   DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
186     It should set the essential boundary data for the local portion of the solution X, as well its time derivative X_t (if it is not NULL).
187     Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.
188 
189   Logically Collective
190 
191   Input Parameters:
192 + dm   - `DM` to associate callback with
193 . func - local function evaluation
194 - ctx  - context for function evaluation
195 
196   Level: intermediate
197 
198   Note:
199   This function is somewhat optional: boundary data could potentially be inserted by a function passed to
200   `DMTSSetIFunctionLocal()`.  The use case for this function is for discretizations with constraints (see
201   `DMGetDefaultConstraints()`): this function inserts boundary values before constraint interpolation.
202 
203 .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
204 @*/
205 PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
206 {
207   DMTS        tdm;
208   DMTS_Local *dmlocalts;
209 
210   PetscFunctionBegin;
211   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
212   PetscCall(DMGetDMTSWrite(dm, &tdm));
213   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
214 
215   dmlocalts->boundarylocal    = func;
216   dmlocalts->boundarylocalctx = ctx;
217 
218   PetscFunctionReturn(PETSC_SUCCESS);
219 }
220 
221 /*@C
222   DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
223       containing the local vector information PLUS ghost point information. It should compute a result for all local
224       elements and `DM` will automatically accumulate the overlapping values.
225 
226   Logically Collective
227 
228   Input Parameter:
229 . dm   - `DM` to associate callback with
230 
231   Output Parameters:
232 + func - local function evaluation
233 - ctx  - context for function evaluation
234 
235   Level: beginner
236 
237 .seealso: [](ch_ts), `DM`, `DMTSSetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
238 @*/
239 PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx)
240 {
241   DMTS           tdm;
242   DMTS_Local    *dmlocalts;
243   PetscErrorCode ierr;
244 
245   PetscFunctionBegin;
246   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
247   ierr = DMGetDMTS(dm, &tdm);
248   CHKERRQ(ierr);
249   ierr = DMLocalTSGetContext(dm, tdm, &dmlocalts);
250   CHKERRQ(ierr);
251   if (func) {
252     PetscValidPointer(func, 2);
253     *func = dmlocalts->ifunctionlocal;
254   }
255   if (ctx) {
256     PetscValidPointer(ctx, 3);
257     *ctx = dmlocalts->ifunctionlocalctx;
258   }
259   PetscFunctionReturn(PETSC_SUCCESS);
260 }
261 
262 /*@C
263   DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
264       containing the local vector information PLUS ghost point information. It should compute a result for all local
265       elements and `DM` will automatically accumulate the overlapping values.
266 
267   Logically Collective
268 
269   Input Parameters:
270 + dm   - `DM` to associate callback with
271 . func - local function evaluation
272 - ctx  - context for function evaluation
273 
274   Level: beginner
275 
276 .seealso: [](ch_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
277 @*/
278 PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
279 {
280   DMTS        tdm;
281   DMTS_Local *dmlocalts;
282 
283   PetscFunctionBegin;
284   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
285   PetscCall(DMGetDMTSWrite(dm, &tdm));
286   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
287 
288   dmlocalts->ifunctionlocal    = func;
289   dmlocalts->ifunctionlocalctx = ctx;
290 
291   PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts));
292   if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
293     PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
294   }
295   PetscFunctionReturn(PETSC_SUCCESS);
296 }
297 
298 /*@C
299   DMTSGetIJacobianLocal - get a local Jacobian evaluation function
300 
301   Logically Collective
302 
303   Input Parameter:
304 . dm - `DM` to associate callback with
305 
306   Output Parameters:
307 + func - local Jacobian evaluation
308 - ctx - optional context for local Jacobian evaluation
309 
310   Level: beginner
311 
312 .seealso: [](ch_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
313 @*/
314 PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void **ctx)
315 {
316   DMTS           tdm;
317   DMTS_Local    *dmlocalts;
318   PetscErrorCode ierr;
319 
320   PetscFunctionBegin;
321   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
322   ierr = DMGetDMTS(dm, &tdm);
323   CHKERRQ(ierr);
324   ierr = DMLocalTSGetContext(dm, tdm, &dmlocalts);
325   CHKERRQ(ierr);
326   if (func) {
327     PetscValidPointer(func, 2);
328     *func = dmlocalts->ijacobianlocal;
329   }
330   if (ctx) {
331     PetscValidPointer(ctx, 3);
332     *ctx = dmlocalts->ijacobianlocalctx;
333   }
334   PetscFunctionReturn(PETSC_SUCCESS);
335 }
336 
337 /*@C
338   DMTSSetIJacobianLocal - set a local Jacobian evaluation function
339 
340   Logically Collective
341 
342   Input Parameters:
343 + dm - `DM` to associate callback with
344 . func - local Jacobian evaluation
345 - ctx - optional context for local Jacobian evaluation
346 
347   Level: beginner
348 
349 .seealso: [](ch_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
350 @*/
351 PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
352 {
353   DMTS        tdm;
354   DMTS_Local *dmlocalts;
355 
356   PetscFunctionBegin;
357   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
358   PetscCall(DMGetDMTSWrite(dm, &tdm));
359   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
360 
361   dmlocalts->ijacobianlocal    = func;
362   dmlocalts->ijacobianlocalctx = ctx;
363 
364   PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
365   PetscFunctionReturn(PETSC_SUCCESS);
366 }
367 
368 /*@C
369   DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
370       containing the local vector information PLUS ghost point information. It should compute a result for all local
371       elements and `DM` will automatically accumulate the overlapping values.
372 
373   Logically Collective
374 
375   Input Parameter:
376 . dm   - `DM` to associate callback with
377 
378   Output Parameters:
379 + func - local function evaluation
380 - ctx  - context for function evaluation
381 
382   Level: beginner
383 
384 .seealso: [](ch_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
385 @*/
386 PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx)
387 {
388   DMTS           tdm;
389   DMTS_Local    *dmlocalts;
390   PetscErrorCode ierr;
391 
392   PetscFunctionBegin;
393   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
394   ierr = DMGetDMTS(dm, &tdm);
395   CHKERRQ(ierr);
396   ierr = DMLocalTSGetContext(dm, tdm, &dmlocalts);
397   CHKERRQ(ierr);
398   if (func) {
399     PetscValidPointer(func, 2);
400     *func = dmlocalts->rhsfunctionlocal;
401   }
402   if (ctx) {
403     PetscValidPointer(ctx, 3);
404     *ctx = dmlocalts->rhsfunctionlocalctx;
405   }
406   PetscFunctionReturn(PETSC_SUCCESS);
407 }
408 
409 /*@C
410   DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
411       containing the local vector information PLUS ghost point information. It should compute a result for all local
412       elements and `DM` will automatically accumulate the overlapping values.
413 
414   Logically Collective
415 
416   Input Parameters:
417 + dm   - `DM` to associate callback with
418 . func - local function evaluation
419 - ctx  - context for function evaluation
420 
421   Level: beginner
422 
423 .seealso: [](ch_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
424 @*/
425 PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
426 {
427   DMTS        tdm;
428   DMTS_Local *dmlocalts;
429 
430   PetscFunctionBegin;
431   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
432   PetscCall(DMGetDMTSWrite(dm, &tdm));
433   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
434 
435   dmlocalts->rhsfunctionlocal    = func;
436   dmlocalts->rhsfunctionlocalctx = ctx;
437 
438   PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts));
439   PetscFunctionReturn(PETSC_SUCCESS);
440 }
441 
442 /*@C
443   DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.
444 
445   Collective
446 
447   Input Parameter:
448 . dm   - `DM` providing the mass matrix
449 
450   Level: developer
451 
452   Note:
453   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.
454 
455 .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
456 @*/
457 PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
458 {
459   DMTS        tdm;
460   DMTS_Local *dmlocalts;
461   const char *prefix;
462 
463   PetscFunctionBegin;
464   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
465   PetscCall(DMGetDMTSWrite(dm, &tdm));
466   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
467   PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass));
468   PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass));
469   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
470   PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix));
471   PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_"));
472   PetscCall(KSPSetFromOptions(dmlocalts->kspmass));
473   PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass));
474   PetscFunctionReturn(PETSC_SUCCESS);
475 }
476 
477 /*@C
478   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.
479 
480   Collective
481 
482   Input Parameter:
483 . dm   - `DM` providing the mass matrix
484 
485   Level: developer
486 
487   Note:
488   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.
489   Since the matrix is lumped, inversion is trivial.
490 
491 .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
492 @*/
493 PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
494 {
495   DMTS        tdm;
496   DMTS_Local *dmlocalts;
497 
498   PetscFunctionBegin;
499   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
500   PetscCall(DMGetDMTSWrite(dm, &tdm));
501   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
502   PetscCall(DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv));
503   PetscCall(VecReciprocal(dmlocalts->lumpedmassinv));
504   PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view"));
505   PetscFunctionReturn(PETSC_SUCCESS);
506 }
507 
508 /*@C
509   DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist.
510 
511   Logically Collective
512 
513   Input Parameter:
514 . dm   - `DM` providing the mass matrix
515 
516   Level: developer
517 
518 .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
519 @*/
520 PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
521 {
522   DMTS        tdm;
523   DMTS_Local *dmlocalts;
524 
525   PetscFunctionBegin;
526   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
527   PetscCall(DMGetDMTSWrite(dm, &tdm));
528   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
529   PetscCall(VecDestroy(&dmlocalts->lumpedmassinv));
530   PetscCall(MatDestroy(&dmlocalts->mass));
531   PetscCall(KSPDestroy(&dmlocalts->kspmass));
532   PetscFunctionReturn(PETSC_SUCCESS);
533 }
534