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) {
114:     PetscCall(VecPointwiseMult(F, dmlocalts->lumpedmassinv, F));
115:   } else if (dmlocalts->kspmass) {
116:     Vec tmp;

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

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

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

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

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

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

194:   Logically Collective

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

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

208:   Level: intermediate

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

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

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

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

229:   PetscFunctionBegin;
231:   PetscCall(DMGetDMTSWrite(dm, &tdm));
232:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

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

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

244:   Logically Collective

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

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

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

261:   Level: beginner

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

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

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

290:   Logically Collective

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

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

305:   Level: beginner

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

314:   PetscFunctionBegin;
316:   PetscCall(DMGetDMTSWrite(dm, &tdm));
317:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

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

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

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

332:   Logically Collective

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

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

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

351:   Level: beginner

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

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

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

378:   Logically Collective

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

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

395:   Level: beginner

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

404:   PetscFunctionBegin;
406:   PetscCall(DMGetDMTSWrite(dm, &tdm));
407:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

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

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

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

421:   Logically Collective

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

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

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

437:   Level: beginner

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

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

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

466:   Logically Collective

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

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

480:   Level: beginner

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

489:   PetscFunctionBegin;
491:   PetscCall(DMGetDMTSWrite(dm, &tdm));
492:   PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));

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

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

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

504:   Collective

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

509:   Level: developer

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

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

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

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

539:   Collective

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

544:   Level: developer

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

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

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

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

570:   Logically Collective

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

575:   Level: developer

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

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