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