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, void *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));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
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;
88: PetscFunctionBegin;
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;
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: }
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;
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;
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));
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: }
184: /*@C
185: DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
187: Logically Collective
189: Input Parameters:
190: + dm - `DM` to associate callback with
191: . func - local function evaluation
192: - ctx - context for function evaluation
194: Level: intermediate
196: Notes:
197: `func` should set the essential boundary data for the local portion of the solution, as
198: well its time derivative (if it is not `NULL`).
200: Vectors are initialized to zero before this function, so it is only needed for non
201: homogeneous data.
203: This function is somewhat optional: boundary data could potentially be inserted by a function
204: passed to `DMTSSetIFunctionLocal()`. The use case for this function is for discretizations
205: with constraints (see `DMGetDefaultConstraints()`): this function inserts boundary values
206: before constraint interpolation.
208: .seealso: [](ch_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
209: @*/
210: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
211: {
212: DMTS tdm;
213: DMTS_Local *dmlocalts;
215: PetscFunctionBegin;
217: PetscCall(DMGetDMTSWrite(dm, &tdm));
218: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
220: dmlocalts->boundarylocal = func;
221: dmlocalts->boundarylocalctx = ctx;
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: /*@C
227: DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
228: containing the local vector information PLUS ghost point information. It should compute a result for all local
229: elements and `DM` will automatically accumulate the overlapping values.
231: Logically Collective
233: Input Parameter:
234: . dm - `DM` to associate callback with
236: Output Parameters:
237: + func - local function evaluation
238: - ctx - context for function evaluation
240: Level: beginner
242: .seealso: [](ch_ts), `DM`, `DMTSSetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
243: @*/
244: PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx)
245: {
246: DMTS tdm;
247: DMTS_Local *dmlocalts;
249: PetscFunctionBegin;
251: PetscCall(DMGetDMTS(dm, &tdm));
252: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
253: if (func) {
254: PetscAssertPointer(func, 2);
255: *func = dmlocalts->ifunctionlocal;
256: }
257: if (ctx) {
258: PetscAssertPointer(ctx, 3);
259: *ctx = dmlocalts->ifunctionlocalctx;
260: }
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@C
265: DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
266: containing the local vector information PLUS ghost point information. It should compute a result for all local
267: elements and `DM` will automatically accumulate the overlapping values.
269: Logically Collective
271: Input Parameters:
272: + dm - `DM` to associate callback with
273: . func - local function evaluation
274: - ctx - context for function evaluation
276: Level: beginner
278: .seealso: [](ch_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
279: @*/
280: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
281: {
282: DMTS tdm;
283: DMTS_Local *dmlocalts;
285: PetscFunctionBegin;
287: PetscCall(DMGetDMTSWrite(dm, &tdm));
288: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
290: dmlocalts->ifunctionlocal = func;
291: dmlocalts->ifunctionlocalctx = ctx;
293: PetscCall(DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts));
294: if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
295: PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
296: }
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: /*@C
301: DMTSGetIJacobianLocal - get a local Jacobian evaluation function
303: Logically Collective
305: Input Parameter:
306: . dm - `DM` to associate callback with
308: Output Parameters:
309: + func - local Jacobian evaluation
310: - ctx - optional context for local Jacobian evaluation
312: Level: beginner
314: .seealso: [](ch_ts), `DM`, `DMTSSetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
315: @*/
316: PetscErrorCode DMTSGetIJacobianLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void **ctx)
317: {
318: DMTS tdm;
319: DMTS_Local *dmlocalts;
321: PetscFunctionBegin;
323: PetscCall(DMGetDMTS(dm, &tdm));
324: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
325: if (func) {
326: PetscAssertPointer(func, 2);
327: *func = dmlocalts->ijacobianlocal;
328: }
329: if (ctx) {
330: PetscAssertPointer(ctx, 3);
331: *ctx = dmlocalts->ijacobianlocalctx;
332: }
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@C
337: DMTSSetIJacobianLocal - set a local Jacobian evaluation function
339: Logically Collective
341: Input Parameters:
342: + dm - `DM` to associate callback with
343: . func - local Jacobian evaluation
344: - ctx - optional context for local Jacobian evaluation
346: Level: beginner
348: .seealso: [](ch_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
349: @*/
350: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
351: {
352: DMTS tdm;
353: DMTS_Local *dmlocalts;
355: PetscFunctionBegin;
357: PetscCall(DMGetDMTSWrite(dm, &tdm));
358: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
360: dmlocalts->ijacobianlocal = func;
361: dmlocalts->ijacobianlocalctx = ctx;
363: PetscCall(DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts));
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@C
368: DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
369: containing the local vector information PLUS ghost point information. It should compute a result for all local
370: elements and `DM` will automatically accumulate the overlapping values.
372: Logically Collective
374: Input Parameter:
375: . dm - `DM` to associate callback with
377: Output Parameters:
378: + func - local function evaluation
379: - ctx - context for function evaluation
381: Level: beginner
383: .seealso: [](ch_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
384: @*/
385: PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx)
386: {
387: DMTS tdm;
388: DMTS_Local *dmlocalts;
390: PetscFunctionBegin;
392: PetscCall(DMGetDMTS(dm, &tdm));
393: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
394: if (func) {
395: PetscAssertPointer(func, 2);
396: *func = dmlocalts->rhsfunctionlocal;
397: }
398: if (ctx) {
399: PetscAssertPointer(ctx, 3);
400: *ctx = dmlocalts->rhsfunctionlocalctx;
401: }
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: /*@C
406: DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
407: containing the local vector information PLUS ghost point information. It should compute a result for all local
408: elements and `DM` will automatically accumulate the overlapping values.
410: Logically Collective
412: Input Parameters:
413: + dm - `DM` to associate callback with
414: . func - local function evaluation
415: - ctx - context for function evaluation
417: Level: beginner
419: .seealso: [](ch_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
420: @*/
421: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
422: {
423: DMTS tdm;
424: DMTS_Local *dmlocalts;
426: PetscFunctionBegin;
428: PetscCall(DMGetDMTSWrite(dm, &tdm));
429: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
431: dmlocalts->rhsfunctionlocal = func;
432: dmlocalts->rhsfunctionlocalctx = ctx;
434: PetscCall(DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@C
439: DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.
441: Collective
443: Input Parameter:
444: . dm - `DM` providing the mass matrix
446: Level: developer
448: Note:
449: 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.
451: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
452: @*/
453: PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
454: {
455: DMTS tdm;
456: DMTS_Local *dmlocalts;
457: const char *prefix;
459: PetscFunctionBegin;
461: PetscCall(DMGetDMTSWrite(dm, &tdm));
462: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
463: PetscCall(DMCreateMassMatrix(dm, dm, &dmlocalts->mass));
464: PetscCall(KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass));
465: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
466: PetscCall(KSPSetOptionsPrefix(dmlocalts->kspmass, prefix));
467: PetscCall(KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_"));
468: PetscCall(KSPSetFromOptions(dmlocalts->kspmass));
469: PetscCall(KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass));
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: /*@C
474: 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.
476: Collective
478: Input Parameter:
479: . dm - `DM` providing the mass matrix
481: Level: developer
483: Note:
484: 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.
485: Since the matrix is lumped, inversion is trivial.
487: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
488: @*/
489: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
490: {
491: DMTS tdm;
492: DMTS_Local *dmlocalts;
494: PetscFunctionBegin;
496: PetscCall(DMGetDMTSWrite(dm, &tdm));
497: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
498: PetscCall(DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv));
499: PetscCall(VecReciprocal(dmlocalts->lumpedmassinv));
500: PetscCall(VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view"));
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: /*@C
505: DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist.
507: Logically Collective
509: Input Parameter:
510: . dm - `DM` providing the mass matrix
512: Level: developer
514: .seealso: [](ch_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMTS`
515: @*/
516: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
517: {
518: DMTS tdm;
519: DMTS_Local *dmlocalts;
521: PetscFunctionBegin;
523: PetscCall(DMGetDMTSWrite(dm, &tdm));
524: PetscCall(DMLocalTSGetContext(dm, tdm, &dmlocalts));
525: PetscCall(VecDestroy(&dmlocalts->lumpedmassinv));
526: PetscCall(MatDestroy(&dmlocalts->mass));
527: PetscCall(KSPDestroy(&dmlocalts->kspmass));
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }