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 (*ifunctionpre)(DM, PetscReal, Vec, Vec, void *);
7: PetscErrorCode (*ifunctionlocal)(DM, PetscReal, Vec, Vec, Vec, void *);
8: PetscErrorCode (*ijacobianlocal)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
9: PetscErrorCode (*rhsfunctionlocal)(DM, PetscReal, Vec, Vec, void *);
10: void *boundarylocalctx;
11: void *ifunctionprectx;
12: void *ifunctionlocalctx;
13: void *ijacobianlocalctx;
14: void *rhsfunctionlocalctx;
15: Vec lumpedmassinv;
16: Mat mass;
17: KSP kspmass;
18: } DMTS_Local;
20: static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm)
21: {
22: PetscFunctionBegin;
23: PetscCall(PetscFree(tdm->data));
24: PetscFunctionReturn(PETSC_SUCCESS);
25: }
27: static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
28: {
29: PetscFunctionBegin;
30: PetscCall(PetscNew((DMTS_Local **)&tdm->data));
31: if (oldtdm->data) PetscCall(PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local)));
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts)
36: {
37: PetscFunctionBegin;
38: *dmlocalts = NULL;
39: if (!tdm->data) {
40: PetscCall(PetscNew((DMTS_Local **)&tdm->data));
42: tdm->ops->destroy = DMTSDestroy_DMLocal;
43: tdm->ops->duplicate = DMTSDuplicate_DMLocal;
44: }
45: *dmlocalts = (DMTS_Local *)tdm->data;
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, PetscCtx ctx)
50: {
51: DM dm;
52: Vec locX, locX_t, locF;
53: DMTS_Local *dmlocalts = (DMTS_Local *)ctx;
55: PetscFunctionBegin;
60: PetscCall(TSGetDM(ts, &dm));
61: PetscCall(DMGetLocalVector(dm, &locX));
62: PetscCall(DMGetLocalVector(dm, &locX_t));
63: PetscCall(DMGetLocalVector(dm, &locF));
64: PetscCall(VecZeroEntries(locX));
65: PetscCall(VecZeroEntries(locX_t));
66: if (dmlocalts->ifunctionpre) PetscCall((*dmlocalts->ifunctionpre)(dm, time, X, X_t, dmlocalts->ifunctionprectx));
67: if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx));
68: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
69: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
70: PetscCall(DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t));
71: PetscCall(DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t));
72: PetscCall(VecZeroEntries(locF));
73: CHKMEMQ;
74: PetscCall((*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx));
75: CHKMEMQ;
76: PetscCall(VecZeroEntries(F));
77: PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F));
78: PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F));
79: PetscCall(DMRestoreLocalVector(dm, &locX));
80: PetscCall(DMRestoreLocalVector(dm, &locX_t));
81: PetscCall(DMRestoreLocalVector(dm, &locF));
83: /* remove nullspace from residual */
84: {
85: MatNullSpace nullsp;
86: PetscCall(PetscObjectQuery((PetscObject)dm, "__dmtsnullspace", (PetscObject *)&nullsp));
87: if (nullsp) PetscCall(MatNullSpaceRemove(nullsp, F));
88: }
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, PetscCtx ctx)
93: {
94: DM dm;
95: Vec locX, locF;
96: DMTS_Local *dmlocalts = (DMTS_Local *)ctx;
98: PetscFunctionBegin;
102: PetscCall(TSGetDM(ts, &dm));
103: PetscCall(DMGetLocalVector(dm, &locX));
104: PetscCall(DMGetLocalVector(dm, &locF));
105: PetscCall(VecZeroEntries(locX));
106: if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, NULL, dmlocalts->boundarylocalctx));
107: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
108: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
109: PetscCall(VecZeroEntries(locF));
110: CHKMEMQ;
111: PetscCall((*dmlocalts->rhsfunctionlocal)(dm, time, locX, locF, dmlocalts->rhsfunctionlocalctx));
112: CHKMEMQ;
113: PetscCall(VecZeroEntries(F));
114: PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F));
115: PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F));
116: if (dmlocalts->lumpedmassinv) PetscCall(VecPointwiseMult(F, dmlocalts->lumpedmassinv, F));
117: else if (dmlocalts->kspmass) {
118: Vec tmp;
120: PetscCall(DMGetGlobalVector(dm, &tmp));
121: PetscCall(KSPSolve(dmlocalts->kspmass, F, tmp));
122: PetscCall(VecCopy(tmp, F));
123: PetscCall(DMRestoreGlobalVector(dm, &tmp));
124: }
125: PetscCall(DMRestoreLocalVector(dm, &locX));
126: PetscCall(DMRestoreLocalVector(dm, &locF));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, PetscCtx ctx)
131: {
132: DM dm;
133: Vec locX, locX_t;
134: DMTS_Local *dmlocalts = (DMTS_Local *)ctx;
136: PetscFunctionBegin;
137: PetscCall(TSGetDM(ts, &dm));
138: if (dmlocalts->ijacobianlocal) {
139: PetscCall(DMGetLocalVector(dm, &locX));
140: PetscCall(DMGetLocalVector(dm, &locX_t));
141: PetscCall(VecZeroEntries(locX));
142: PetscCall(VecZeroEntries(locX_t));
143: if (dmlocalts->boundarylocal) PetscCall((*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx));
144: PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX));
145: PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX));
146: PetscCall(DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t));
147: PetscCall(DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t));
148: CHKMEMQ;
149: PetscCall((*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx));
150: CHKMEMQ;
151: PetscCall(DMRestoreLocalVector(dm, &locX));
152: PetscCall(DMRestoreLocalVector(dm, &locX_t));
153: } else {
154: MatFDColoring fdcoloring;
155: PetscCall(PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring));
156: if (!fdcoloring) {
157: ISColoring coloring;
159: PetscCall(DMCreateColoring(dm, dm->coloringtype, &coloring));
160: PetscCall(MatFDColoringCreate(B, coloring, &fdcoloring));
161: PetscCall(ISColoringDestroy(&coloring));
162: switch (dm->coloringtype) {
163: case IS_COLORING_GLOBAL:
164: PetscCall(MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn *)(PetscVoidFn *)TSComputeIFunction_DMLocal, dmlocalts));
165: break;
166: default:
167: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
168: }
169: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix));
170: PetscCall(MatFDColoringSetFromOptions(fdcoloring));
171: PetscCall(MatFDColoringSetUp(B, coloring, fdcoloring));
172: PetscCall(PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring));
173: PetscCall(PetscObjectDereference((PetscObject)fdcoloring));
175: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
176: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
177: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
178: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
179: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
180: */
181: PetscCall(PetscObjectDereference((PetscObject)dm));
182: }
183: PetscCall(MatFDColoringApply(B, fdcoloring, X, ts));
184: }
185: /* This will be redundant if the user called both, but it's too common to forget. */
186: if (A != B) {
187: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
188: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /*@C
194: DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
196: Logically Collective
198: Input Parameters:
199: + dm - `DM` to associate callback with
200: . func - local function evaluation
201: - ctx - context for function evaluation
203: Calling sequence of `func`:
204: + dm - the `DM`
205: . t - the current time
206: . u - the current solution
207: . f - output, the computed right hand side function
208: - ctx - the application context for the function
210: Level: intermediate
212: Notes:
213: `func` should set the essential boundary data for the local portion of the solution, as
214: well its time derivative (if it is not `NULL`).
216: Vectors are initialized to zero before this function, so it is only needed for non
217: homogeneous data.
219: This function is somewhat optional: boundary data could potentially be inserted by a function
220: passed to `DMTSSetIFunctionLocal()`. The use case for this function is for discretizations
221: with constraints (see `DMGetDefaultConstraints()`): this function inserts boundary values
222: before constraint interpolation.
224: .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
225: @*/
226: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec f, PetscCtx ctx), PetscCtx ctx)
227: {
228: DMTS tdm;
229: DMTS_Local *dmlocalts;
231: PetscFunctionBegin;
233: PetscCall(DMGetDMTSWrite(dm, &tdm));
234: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
236: dmlocalts->boundarylocal = func;
237: dmlocalts->boundarylocalctx = ctx;
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /*@C
242: DMTSSetIFunctionPre - set a pre-evaluation callback for each local implicit function evaluation. The callback function provided is called at the beginning of `TSComputeIFunction()` before the function provided with `TSSetIFunctionLocal()` and `DMTSSetBoundaryLocal()` is called.
244: Logically Collective
246: Input Parameters:
247: + dm - `DM` to associate callback with
248: . func - function evaluation
249: - ctx - context for function evaluation
251: Calling sequence of `func`:
252: + dm - the `DM`
253: . time - the current time
254: . u - the global solution
255: . u_t - the global time derivative
256: - ctx - the user context
258: Level: intermediate
260: Notes:
261: `func` should perform any setup required before the local implicit function evaluation provided by `TSSetIFunctionLocal()` can be performed, such as computing auxiliary data via a subsolve. This function is necessary because some setups may require global communication, such as the Poisson solve in the Vlasov-Poisson system, and cannot be embedded in only a local evaluation, yet we want to specify the local callback for the main equations using `TSSetIFunctionLocal()`.
263: .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
264: @*/
265: PetscErrorCode DMTSSetIFunctionPre(DM dm, PetscErrorCode (*func)(DM dm, PetscReal time, Vec u, Vec u_t, void *ctx), void *ctx)
266: {
267: DMTS tdm;
268: DMTS_Local *dmlocalts;
270: PetscFunctionBegin;
272: PetscCall(DMGetDMTSWrite(dm, &tdm));
273: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
275: dmlocalts->ifunctionpre = func;
276: dmlocalts->ifunctionprectx = ctx;
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /*@C
281: DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
282: containing the local vector information PLUS ghost point information. It should compute a result for all local
283: elements and `DM` will automatically accumulate the overlapping values.
285: Logically Collective
287: Input Parameter:
288: . dm - `DM` to associate callback with
290: Output Parameters:
291: + func - local function evaluation
292: - ctx - context for function evaluation
294: Calling sequence of `func`:
295: + dm - the `DM`
296: . t - the current time
297: . u - the current solution
298: . udot - the derivative of `u`
299: . F - output, the computed implicit function
300: - ctx - the application context for the function
302: Level: beginner
304: .seealso: [](ch_ts), `DM`, `DMTSSetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
305: @*/
306: PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM dm, PetscReal t, Vec u, Vec udot, Vec F, PetscCtx ctx), PetscCtxRt ctx)
307: {
308: DMTS tdm;
309: DMTS_Local *dmlocalts;
311: PetscFunctionBegin;
313: PetscCall(DMGetDMTS(dm, &tdm));
314: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
315: if (func) {
316: PetscAssertPointer(func, 2);
317: *func = dmlocalts->ifunctionlocal;
318: }
319: if (ctx) {
320: PetscAssertPointer(ctx, 3);
321: *(void **)ctx = dmlocalts->ifunctionlocalctx;
322: }
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /*@C
327: DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
328: containing the local vector information PLUS ghost point information. It should compute a result for all local
329: elements and `DM` will automatically accumulate the overlapping values.
331: Logically Collective
333: Input Parameters:
334: + dm - `DM` to associate callback with
335: . func - local function evaluation
336: - ctx - context for function evaluation
338: Calling sequence of `func`:
339: + dm - the `DM`
340: . t - the current time
341: . u - the current solution
342: . udot - the derivative of `u`
343: . F - output, the computed implicit function
344: - ctx - the application context for the function
346: Level: beginner
348: .seealso: [](ch_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
349: @*/
350: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec udot, Vec F, PetscCtx ctx), PetscCtx ctx)
351: {
352: DMTS tdm;
353: DMTS_Local *dmlocalts;
355: PetscFunctionBegin;
357: PetscCall(DMGetDMTSWrite(dm, &tdm));
358: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
360: dmlocalts->ifunctionlocal = func;
361: dmlocalts->ifunctionlocalctx = ctx;
363: PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts));
364: if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
365: PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
366: }
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: /*@C
371: DMTSGetIJacobianLocal - get a local Jacobian evaluation function
373: Logically Collective
375: Input Parameter:
376: . dm - `DM` to associate callback with
378: Output Parameters:
379: + func - local Jacobian evaluation
380: - ctx - optional context for local Jacobian evaluation
382: Calling sequence of `func`:
383: + dm - the `DM`
384: . t - the current time
385: . u - the current solution
386: . udot - the derivative of `u`
387: . shift - the shift factoring arising from the implicit time-step
388: . J - output, the Jacobian
389: . Jpre - output, matrix from which to compute the preconditioner for `J`, often the same as `J`
390: - ctx - the application context for the function
392: Level: beginner
394: .seealso: [](ch_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
395: @*/
396: PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM dm, PetscReal t, Vec u, Vec udot, PetscReal shift, Mat J, Mat Jpre, PetscCtx ctx), PetscCtxRt ctx)
397: {
398: DMTS tdm;
399: DMTS_Local *dmlocalts;
401: PetscFunctionBegin;
403: PetscCall(DMGetDMTS(dm, &tdm));
404: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
405: if (func) {
406: PetscAssertPointer(func, 2);
407: *func = dmlocalts->ijacobianlocal;
408: }
409: if (ctx) {
410: PetscAssertPointer(ctx, 3);
411: *(void **)ctx = dmlocalts->ijacobianlocalctx;
412: }
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: /*@C
417: DMTSSetIJacobianLocal - set a local Jacobian evaluation function
419: Logically Collective
421: Input Parameters:
422: + dm - `DM` to associate callback with
423: . func - local Jacobian evaluation
424: - ctx - optional context for local Jacobian evaluation
426: Calling sequence of `func`:
427: + dm - the `DM`
428: . t - the current time
429: . u - the current solution
430: . udot - the derivative of `u`
431: . shift - the shift factoring arising from the implicit time-step
432: . J - output, the Jacobian
433: . Jpre - output, matrix from which to compute the preconditioner for `J`, often the same as `J`
434: - ctx - the application context for the function
436: Level: beginner
438: .seealso: [](ch_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
439: @*/
440: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec udot, PetscReal shift, Mat J, Mat Jpre, PetscCtx ctx), PetscCtx ctx)
441: {
442: DMTS tdm;
443: DMTS_Local *dmlocalts;
445: PetscFunctionBegin;
447: PetscCall(DMGetDMTSWrite(dm, &tdm));
448: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
450: dmlocalts->ijacobianlocal = func;
451: dmlocalts->ijacobianlocalctx = ctx;
453: PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@C
458: DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
459: containing the local vector information PLUS ghost point information. It should compute a result for all local
460: elements and `DM` will automatically accumulate the overlapping values.
462: Logically Collective
464: Input Parameter:
465: . dm - `DM` to associate callback with
467: Output Parameters:
468: + func - local function evaluation
469: - ctx - context for function evaluation
471: Calling sequence of `func`:
472: + dm - the `DM`
473: . t - the current time
474: . u - the current solution
475: . udot - output, the evaluated right hand side
476: - ctx - the application context for the function
478: Level: beginner
480: .seealso: [](ch_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
481: @*/
482: PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM dm, PetscReal t, Vec u, Vec udot, PetscCtx ctx), PetscCtxRt ctx)
483: {
484: DMTS tdm;
485: DMTS_Local *dmlocalts;
487: PetscFunctionBegin;
489: PetscCall(DMGetDMTS(dm, &tdm));
490: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
491: if (func) {
492: PetscAssertPointer(func, 2);
493: *func = dmlocalts->rhsfunctionlocal;
494: }
495: if (ctx) {
496: PetscAssertPointer(ctx, 3);
497: *(void **)ctx = dmlocalts->rhsfunctionlocalctx;
498: }
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /*@C
503: DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
504: containing the local vector information PLUS ghost point information. It should compute a result for all local
505: elements and `DM` will automatically accumulate the overlapping values.
507: Logically Collective
509: Input Parameters:
510: + dm - `DM` to associate callback with
511: . func - local function evaluation
512: - ctx - context for function evaluation
514: Calling sequence of `func`:
515: + dm - the `DM`
516: . t - the current time
517: . u - the current solution
518: . f - output, the evaluated right hand side
519: - ctx - the application context for the function
521: Level: beginner
523: .seealso: [](ch_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
524: @*/
525: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM dm, PetscReal t, Vec u, Vec f, PetscCtx ctx), PetscCtx ctx)
526: {
527: DMTS tdm;
528: DMTS_Local *dmlocalts;
530: PetscFunctionBegin;
532: PetscCall(DMGetDMTSWrite(dm, &tdm));
533: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
535: dmlocalts->rhsfunctionlocal = func;
536: dmlocalts->rhsfunctionlocalctx = ctx;
538: PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*@
543: DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.
545: Collective
547: Input Parameter:
548: . dm - `DM` providing the mass matrix
550: Level: developer
552: Note:
553: 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.
555: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
556: @*/
557: PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
558: {
559: DMTS tdm;
560: DMTS_Local *dmlocalts;
561: const char *prefix;
563: PetscFunctionBegin;
565: PetscCall(DMGetDMTSWrite(dm, &tdm));
566: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
567: PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass));
568: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass));
569: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
570: PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix));
571: PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_"));
572: PetscCall(KSPSetFromOptions(dmlocalts->kspmass));
573: PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass));
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@
578: 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.
580: Collective
582: Input Parameter:
583: . dm - `DM` providing the mass matrix
585: Level: developer
587: Note:
588: 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.
589: Since the matrix is lumped, inversion is trivial.
591: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
592: @*/
593: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
594: {
595: DMTS tdm;
596: DMTS_Local *dmlocalts;
598: PetscFunctionBegin;
600: PetscCall(DMGetDMTSWrite(dm, &tdm));
601: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
602: PetscCall(DMCreateMassMatrixLumped(dm, NULL, &dmlocalts->lumpedmassinv));
603: PetscCall(VecReciprocal(dmlocalts->lumpedmassinv));
604: PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view"));
605: PetscFunctionReturn(PETSC_SUCCESS);
606: }
608: /*@
609: DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist.
611: Logically Collective
613: Input Parameter:
614: . dm - `DM` providing the mass matrix
616: Level: developer
618: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMTS`
619: @*/
620: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
621: {
622: DMTS tdm;
623: DMTS_Local *dmlocalts;
625: PetscFunctionBegin;
627: PetscCall(DMGetDMTSWrite(dm, &tdm));
628: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
629: PetscCall(VecDestroy(&dmlocalts->lumpedmassinv));
630: PetscCall(MatDestroy(&dmlocalts->mass));
631: PetscCall(KSPDestroy(&dmlocalts->kspmass));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }