xref: /petsc/src/ts/utils/dmlocalts.c (revision f97672e55eacc8688507b9471cd7ec2664d7f203)
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(0);
23 }
24 
25 static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
26 {
27   PetscFunctionBegin;
28   PetscCall(PetscNewLog(tdm, (DMTS_Local **) &tdm->data));
29   if (oldtdm->data) PetscCall(PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local)));
30   PetscFunctionReturn(0);
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(PetscNewLog(dm, (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(0);
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(0);
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(0);
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: SETERRQ(PetscObjectComm((PetscObject) ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
158       }
159       PetscCall(PetscObjectSetOptionsPrefix((PetscObject) fdcoloring, ((PetscObject) dm)->prefix));
160       PetscCall(MatFDColoringSetFromOptions(fdcoloring));
161       PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
162       PetscCall(PetscObjectCompose((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject) fdcoloring));
163       PetscCall(PetscObjectDereference((PetscObject) fdcoloring));
164 
165       /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
166        * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
167        * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
168        * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
169        * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
170        */
171       PetscCall(PetscObjectDereference((PetscObject) dm));
172     }
173     PetscCall(MatFDColoringApply(B, fdcoloring, X, ts));
174   }
175   /* This will be redundant if the user called both, but it's too common to forget. */
176   if (A != B) {
177     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
178     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
179   }
180   PetscFunctionReturn(0);
181 }
182 
183 /*@C
184   DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
185     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).
186     Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.
187 
188   Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to
189   DMTSSetIFunctionLocal().  The use case for this function is for discretizations with constraints (see
190   DMGetDefaultConstraints()): this function inserts boundary values before constraint interpolation.
191 
192   Logically Collective
193 
194   Input Parameters:
195 + dm   - DM to associate callback with
196 . func - local function evaluation
197 - ctx  - context for function evaluation
198 
199   Level: intermediate
200 
201 .seealso: `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
202 @*/
203 PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
204 {
205   DMTS           tdm;
206   DMTS_Local    *dmlocalts;
207 
208   PetscFunctionBegin;
209   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
210   PetscCall(DMGetDMTSWrite(dm, &tdm));
211   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
212 
213   dmlocalts->boundarylocal    = func;
214   dmlocalts->boundarylocalctx = ctx;
215 
216   PetscFunctionReturn(0);
217 }
218 
219 /*@C
220   DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
221       containing the local vector information PLUS ghost point information. It should compute a result for all local
222       elements and DMTS will automatically accumulate the overlapping values.
223 
224   Logically Collective
225 
226   Input Parameter:
227 . dm   - DM to associate callback with
228 
229   Output Parameters:
230 + func - local function evaluation
231 - ctx  - context for function evaluation
232 
233   Level: beginner
234 
235 .seealso: `DMTSSetIFunctionLocal(()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
236 @*/
237 PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx)
238 {
239   DMTS           tdm;
240   DMTS_Local    *dmlocalts;
241   PetscErrorCode ierr;
242 
243   PetscFunctionBegin;
244   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
245   ierr = DMGetDMTS(dm, &tdm);CHKERRQ(ierr);
246   ierr = DMLocalTSGetContext(dm, tdm, &dmlocalts);CHKERRQ(ierr);
247   if (func) {PetscValidPointer(func, 2); *func = dmlocalts->ifunctionlocal;}
248   if (ctx)  {PetscValidPointer(ctx, 3);  *ctx  = dmlocalts->ifunctionlocalctx;}
249   PetscFunctionReturn(0);
250 }
251 
252 /*@C
253   DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
254       containing the local vector information PLUS ghost point information. It should compute a result for all local
255       elements and DMTS will automatically accumulate the overlapping values.
256 
257   Logically Collective
258 
259   Input Parameters:
260 + dm   - DM to associate callback with
261 . func - local function evaluation
262 - ctx  - context for function evaluation
263 
264   Level: beginner
265 
266 .seealso: `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
267 @*/
268 PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
269 {
270   DMTS           tdm;
271   DMTS_Local    *dmlocalts;
272 
273   PetscFunctionBegin;
274   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
275   PetscCall(DMGetDMTSWrite(dm, &tdm));
276   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
277 
278   dmlocalts->ifunctionlocal    = func;
279   dmlocalts->ifunctionlocalctx = ctx;
280 
281   PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts));
282   if (!tdm->ops->ijacobian) {  /* Call us for the Jacobian too, can be overridden by the user. */
283     PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
284   }
285   PetscFunctionReturn(0);
286 }
287 
288 /*@C
289   DMTSGetIJacobianLocal - get a local Jacobian evaluation function
290 
291   Logically Collective
292 
293   Input Parameter:
294 . dm - DM to associate callback with
295 
296   Output Parameters:
297 + func - local Jacobian evaluation
298 - ctx - optional context for local Jacobian evaluation
299 
300   Level: beginner
301 
302 .seealso: `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
303 @*/
304 PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void **ctx)
305 {
306   DMTS           tdm;
307   DMTS_Local    *dmlocalts;
308   PetscErrorCode ierr;
309 
310   PetscFunctionBegin;
311   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
312   ierr = DMGetDMTS(dm, &tdm);CHKERRQ(ierr);
313   ierr = DMLocalTSGetContext(dm, tdm, &dmlocalts);CHKERRQ(ierr);
314   if (func) {PetscValidPointer(func, 2); *func = dmlocalts->ijacobianlocal;}
315   if (ctx)  {PetscValidPointer(ctx, 3);  *ctx  = dmlocalts->ijacobianlocalctx;}
316   PetscFunctionReturn(0);
317 }
318 
319 /*@C
320   DMTSSetIJacobianLocal - set a local Jacobian evaluation function
321 
322   Logically Collective
323 
324   Input Parameters:
325 + dm - DM to associate callback with
326 . func - local Jacobian evaluation
327 - ctx - optional context for local Jacobian evaluation
328 
329   Level: beginner
330 
331 .seealso: `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
332 @*/
333 PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
334 {
335   DMTS           tdm;
336   DMTS_Local    *dmlocalts;
337 
338   PetscFunctionBegin;
339   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
340   PetscCall(DMGetDMTSWrite(dm, &tdm));
341   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
342 
343   dmlocalts->ijacobianlocal    = func;
344   dmlocalts->ijacobianlocalctx = ctx;
345 
346   PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
347   PetscFunctionReturn(0);
348 }
349 
350 /*@C
351   DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
352       containing the local vector information PLUS ghost point information. It should compute a result for all local
353       elements and DMTS will automatically accumulate the overlapping values.
354 
355   Logically Collective
356 
357   Input Parameter:
358 . dm   - DM to associate callback with
359 
360   Output Parameters:
361 + func - local function evaluation
362 - ctx  - context for function evaluation
363 
364   Level: beginner
365 
366 .seealso: `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
367 @*/
368 PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx)
369 {
370   DMTS           tdm;
371   DMTS_Local    *dmlocalts;
372   PetscErrorCode ierr;
373 
374   PetscFunctionBegin;
375   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
376   ierr = DMGetDMTS(dm, &tdm);CHKERRQ(ierr);
377   ierr = DMLocalTSGetContext(dm, tdm, &dmlocalts);CHKERRQ(ierr);
378   if (func) {PetscValidPointer(func, 2); *func = dmlocalts->rhsfunctionlocal;}
379   if (ctx)  {PetscValidPointer(ctx, 3);  *ctx  = dmlocalts->rhsfunctionlocalctx;}
380   PetscFunctionReturn(0);
381 }
382 
383 /*@C
384   DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
385       containing the local vector information PLUS ghost point information. It should compute a result for all local
386       elements and DMTS will automatically accumulate the overlapping values.
387 
388   Logically Collective
389 
390   Input Parameters:
391 + dm   - DM to associate callback with
392 . func - local function evaluation
393 - ctx  - context for function evaluation
394 
395   Level: beginner
396 
397 .seealso: `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
398 @*/
399 PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
400 {
401   DMTS           tdm;
402   DMTS_Local    *dmlocalts;
403 
404   PetscFunctionBegin;
405   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
406   PetscCall(DMGetDMTSWrite(dm, &tdm));
407   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
408 
409   dmlocalts->rhsfunctionlocal    = func;
410   dmlocalts->rhsfunctionlocalctx = ctx;
411 
412   PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts));
413   PetscFunctionReturn(0);
414 }
415 
416 /*@C
417   DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given DM, and a solver to invert it, and stores them in the DMTS context.
418 
419   Collective on dm
420 
421   Input Parameters:
422 . dm   - DM providing the mass matrix
423 
424   Note: 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.
425 
426   Level: developer
427 
428 .seealso: `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
429 @*/
430 PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
431 {
432   DMTS           tdm;
433   DMTS_Local    *dmlocalts;
434   const char    *prefix;
435 
436   PetscFunctionBegin;
437   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
438   PetscCall(DMGetDMTSWrite(dm, &tdm));
439   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
440   PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass));
441   PetscCall(KSPCreate(PetscObjectComm((PetscObject) dm), &dmlocalts->kspmass));
442   PetscCall(PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix));
443   PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix));
444   PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_"));
445   PetscCall(KSPSetFromOptions(dmlocalts->kspmass));
446   PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass));
447   PetscFunctionReturn(0);
448 }
449 
450 /*@C
451   DMTSCreateRHSMassMatrixLumped - This creates the lumped mass matrix associated with the given DM, and a solver to invert it, and stores them in the DMTS context.
452 
453   Collective on dm
454 
455   Input Parameters:
456 . dm   - DM providing the mass matrix
457 
458   Note: 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.
459   Since the matrix is lumped, inversion is trivial.
460 
461   Level: developer
462 
463 .seealso: `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
464 @*/
465 PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
466 {
467   DMTS           tdm;
468   DMTS_Local    *dmlocalts;
469 
470   PetscFunctionBegin;
471   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
472   PetscCall(DMGetDMTSWrite(dm, &tdm));
473   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
474   PetscCall(DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv));
475   PetscCall(VecReciprocal(dmlocalts->lumpedmassinv));
476   PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view"));
477   PetscFunctionReturn(0);
478 }
479 
480 /*@C
481   DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the DMTS context, if they exist.
482 
483   Logically Collective
484 
485   Input Parameters:
486 . dm   - DM providing the mass matrix
487 
488   Level: developer
489 
490 .seealso: `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
491 @*/
492 PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
493 {
494   DMTS           tdm;
495   DMTS_Local    *dmlocalts;
496 
497   PetscFunctionBegin;
498   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
499   PetscCall(DMGetDMTSWrite(dm, &tdm));
500   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
501   PetscCall(VecDestroy(&dmlocalts->lumpedmassinv));
502   PetscCall(MatDestroy(&dmlocalts->mass));
503   PetscCall(KSPDestroy(&dmlocalts->kspmass));
504   PetscFunctionReturn(0);
505 }
506