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: PetscFree(tdm->data);
21: return 0;
22: }
24: static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
25: {
26: PetscNew((DMTS_Local **)&tdm->data);
27: if (oldtdm->data) PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local));
28: return 0;
29: }
31: static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts)
32: {
33: *dmlocalts = NULL;
34: if (!tdm->data) {
35: PetscNew((DMTS_Local **)&tdm->data);
37: tdm->ops->destroy = DMTSDestroy_DMLocal;
38: tdm->ops->duplicate = DMTSDuplicate_DMLocal;
39: }
40: *dmlocalts = (DMTS_Local *)tdm->data;
41: return 0;
42: }
44: static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, void *ctx)
45: {
46: DM dm;
47: Vec locX, locX_t, locF;
48: DMTS_Local *dmlocalts = (DMTS_Local *)ctx;
54: TSGetDM(ts, &dm);
55: DMGetLocalVector(dm, &locX);
56: DMGetLocalVector(dm, &locX_t);
57: DMGetLocalVector(dm, &locF);
58: VecZeroEntries(locX);
59: VecZeroEntries(locX_t);
60: if (dmlocalts->boundarylocal) (*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx);
61: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
62: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
63: DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
64: DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
65: VecZeroEntries(locF);
66: CHKMEMQ;
67: (*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx);
68: CHKMEMQ;
69: VecZeroEntries(F);
70: DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F);
71: DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F);
72: DMRestoreLocalVector(dm, &locX);
73: DMRestoreLocalVector(dm, &locX_t);
74: DMRestoreLocalVector(dm, &locF);
75: return 0;
76: }
78: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, void *ctx)
79: {
80: DM dm;
81: Vec locX, locF;
82: DMTS_Local *dmlocalts = (DMTS_Local *)ctx;
87: TSGetDM(ts, &dm);
88: DMGetLocalVector(dm, &locX);
89: DMGetLocalVector(dm, &locF);
90: VecZeroEntries(locX);
91: if (dmlocalts->boundarylocal) (*dmlocalts->boundarylocal)(dm, time, locX, NULL, dmlocalts->boundarylocalctx);
92: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
93: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
94: VecZeroEntries(locF);
95: CHKMEMQ;
96: (*dmlocalts->rhsfunctionlocal)(dm, time, locX, locF, dmlocalts->rhsfunctionlocalctx);
97: CHKMEMQ;
98: VecZeroEntries(F);
99: DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F);
100: DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F);
101: if (dmlocalts->lumpedmassinv) {
102: VecPointwiseMult(F, dmlocalts->lumpedmassinv, F);
103: } else if (dmlocalts->kspmass) {
104: Vec tmp;
106: DMGetGlobalVector(dm, &tmp);
107: KSPSolve(dmlocalts->kspmass, F, tmp);
108: VecCopy(tmp, F);
109: DMRestoreGlobalVector(dm, &tmp);
110: }
111: DMRestoreLocalVector(dm, &locX);
112: DMRestoreLocalVector(dm, &locF);
113: return 0;
114: }
116: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, void *ctx)
117: {
118: DM dm;
119: Vec locX, locX_t;
120: DMTS_Local *dmlocalts = (DMTS_Local *)ctx;
122: TSGetDM(ts, &dm);
123: if (dmlocalts->ijacobianlocal) {
124: DMGetLocalVector(dm, &locX);
125: DMGetLocalVector(dm, &locX_t);
126: VecZeroEntries(locX);
127: VecZeroEntries(locX_t);
128: if (dmlocalts->boundarylocal) (*dmlocalts->boundarylocal)(dm, time, locX, locX_t, dmlocalts->boundarylocalctx);
129: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
130: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
131: DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
132: DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
133: CHKMEMQ;
134: (*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx);
135: CHKMEMQ;
136: DMRestoreLocalVector(dm, &locX);
137: DMRestoreLocalVector(dm, &locX_t);
138: } else {
139: MatFDColoring fdcoloring;
140: PetscObjectQuery((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject *)&fdcoloring);
141: if (!fdcoloring) {
142: ISColoring coloring;
144: DMCreateColoring(dm, dm->coloringtype, &coloring);
145: MatFDColoringCreate(B, coloring, &fdcoloring);
146: ISColoringDestroy(&coloring);
147: switch (dm->coloringtype) {
148: case IS_COLORING_GLOBAL:
149: MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))TSComputeIFunction_DMLocal, dmlocalts);
150: break;
151: default:
152: SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
153: }
154: PetscObjectSetOptionsPrefix((PetscObject)fdcoloring, ((PetscObject)dm)->prefix);
155: MatFDColoringSetFromOptions(fdcoloring);
156: MatFDColoringSetUp(B, coloring, fdcoloring);
157: PetscObjectCompose((PetscObject)dm, "DMDASNES_FDCOLORING", (PetscObject)fdcoloring);
158: PetscObjectDereference((PetscObject)fdcoloring);
160: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
161: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
162: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
163: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
164: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
165: */
166: PetscObjectDereference((PetscObject)dm);
167: }
168: MatFDColoringApply(B, fdcoloring, X, ts);
169: }
170: /* This will be redundant if the user called both, but it's too common to forget. */
171: if (A != B) {
172: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
173: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
174: }
175: return 0;
176: }
178: /*@C
179: DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
180: 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).
181: Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.
183: Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to
184: `DMTSSetIFunctionLocal()`. The use case for this function is for discretizations with constraints (see
185: `DMGetDefaultConstraints()`): this function inserts boundary values before constraint interpolation.
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: .seealso: [](chapter_ts), `DM`, `TS`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
197: @*/
198: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
199: {
200: DMTS tdm;
201: DMTS_Local *dmlocalts;
204: DMGetDMTSWrite(dm, &tdm);
205: DMLocalTSGetContext(dm, tdm, &dmlocalts);
207: dmlocalts->boundarylocal = func;
208: dmlocalts->boundarylocalctx = ctx;
210: return 0;
211: }
213: /*@C
214: DMTSGetIFunctionLocal - get the local implicit function evaluation function. This function is called with local vector
215: containing the local vector information PLUS ghost point information. It should compute a result for all local
216: elements and `DM` will automatically accumulate the overlapping values.
218: Logically Collective
220: Input Parameter:
221: . dm - `DM` to associate callback with
223: Output Parameters:
224: + func - local function evaluation
225: - ctx - context for function evaluation
227: Level: beginner
229: .seealso: [](chapter_ts), `DM`, `DMTSSetIFunctionLocal(()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
230: @*/
231: PetscErrorCode DMTSGetIFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, Vec, void *), void **ctx)
232: {
233: DMTS tdm;
234: DMTS_Local *dmlocalts;
238: DMGetDMTS(dm, &tdm);
239: CHKERRQ(ierr);
240: DMLocalTSGetContext(dm, tdm, &dmlocalts);
241: CHKERRQ(ierr);
242: if (func) {
244: *func = dmlocalts->ifunctionlocal;
245: }
246: if (ctx) {
248: *ctx = dmlocalts->ifunctionlocalctx;
249: }
250: return 0;
251: }
253: /*@C
254: DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
255: containing the local vector information PLUS ghost point information. It should compute a result for all local
256: elements and `DM` will automatically accumulate the overlapping values.
258: Logically Collective
260: Input Parameters:
261: + dm - `DM` to associate callback with
262: . func - local function evaluation
263: - ctx - context for function evaluation
265: Level: beginner
267: .seealso: [](chapter_ts), `DM`, `DMTSGetIFunctionLocal()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
268: @*/
269: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
270: {
271: DMTS tdm;
272: DMTS_Local *dmlocalts;
275: DMGetDMTSWrite(dm, &tdm);
276: DMLocalTSGetContext(dm, tdm, &dmlocalts);
278: dmlocalts->ifunctionlocal = func;
279: dmlocalts->ifunctionlocalctx = ctx;
281: DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts);
282: if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
283: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
284: }
285: return 0;
286: }
288: /*@C
289: DMTSGetIJacobianLocal - get a local Jacobian evaluation function
291: Logically Collective
293: Input Parameter:
294: . dm - `DM` to associate callback with
296: Output Parameters:
297: + func - local Jacobian evaluation
298: - ctx - optional context for local Jacobian evaluation
300: Level: beginner
302: .seealso: [](chapter_ts), `DM`, `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;
311: DMGetDMTS(dm, &tdm);
312: CHKERRQ(ierr);
313: DMLocalTSGetContext(dm, tdm, &dmlocalts);
314: CHKERRQ(ierr);
315: if (func) {
317: *func = dmlocalts->ijacobianlocal;
318: }
319: if (ctx) {
321: *ctx = dmlocalts->ijacobianlocalctx;
322: }
323: return 0;
324: }
326: /*@C
327: DMTSSetIJacobianLocal - set a local Jacobian evaluation function
329: Logically Collective
331: Input Parameters:
332: + dm - `DM` to associate callback with
333: . func - local Jacobian evaluation
334: - ctx - optional context for local Jacobian evaluation
336: Level: beginner
338: .seealso: [](chapter_ts), `DM`, `DMTSGetIJacobianLocal()`, `DMTSSetIFunctionLocal()`, `DMTSSetIJacobian()`, `DMTSSetIFunction()`
339: @*/
340: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
341: {
342: DMTS tdm;
343: DMTS_Local *dmlocalts;
346: DMGetDMTSWrite(dm, &tdm);
347: DMLocalTSGetContext(dm, tdm, &dmlocalts);
349: dmlocalts->ijacobianlocal = func;
350: dmlocalts->ijacobianlocalctx = ctx;
352: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
353: return 0;
354: }
356: /*@C
357: DMTSGetRHSFunctionLocal - get a local rhs function evaluation function. This function is called with local vector
358: containing the local vector information PLUS ghost point information. It should compute a result for all local
359: elements and `DM` will automatically accumulate the overlapping values.
361: Logically Collective
363: Input Parameter:
364: . dm - `DM` to associate callback with
366: Output Parameters:
367: + func - local function evaluation
368: - ctx - context for function evaluation
370: Level: beginner
372: .seealso: [](chapter_ts), `DM`, `DMTSSetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
373: @*/
374: PetscErrorCode DMTSGetRHSFunctionLocal(DM dm, PetscErrorCode (**func)(DM, PetscReal, Vec, Vec, void *), void **ctx)
375: {
376: DMTS tdm;
377: DMTS_Local *dmlocalts;
381: DMGetDMTS(dm, &tdm);
382: CHKERRQ(ierr);
383: DMLocalTSGetContext(dm, tdm, &dmlocalts);
384: CHKERRQ(ierr);
385: if (func) {
387: *func = dmlocalts->rhsfunctionlocal;
388: }
389: if (ctx) {
391: *ctx = dmlocalts->rhsfunctionlocalctx;
392: }
393: return 0;
394: }
396: /*@C
397: DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
398: containing the local vector information PLUS ghost point information. It should compute a result for all local
399: elements and `DM` will automatically accumulate the overlapping values.
401: Logically Collective
403: Input Parameters:
404: + dm - `DM` to associate callback with
405: . func - local function evaluation
406: - ctx - context for function evaluation
408: Level: beginner
410: .seealso: [](chapter_ts), `DM`, `DMTSGetRHSFunctionLocal()`, `DMTSSetRHSFunction()`, `DMTSSetIFunction()`, `DMTSSetIJacobianLocal()`
411: @*/
412: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
413: {
414: DMTS tdm;
415: DMTS_Local *dmlocalts;
418: DMGetDMTSWrite(dm, &tdm);
419: DMLocalTSGetContext(dm, tdm, &dmlocalts);
421: dmlocalts->rhsfunctionlocal = func;
422: dmlocalts->rhsfunctionlocalctx = ctx;
424: DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts);
425: return 0;
426: }
428: /*@C
429: DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given `DM`, and a solver to invert it, and stores them in the `DM` context.
431: Collective
433: Input Parameters:
434: . dm - `DM` providing the mass matrix
436: Level: developer
438: Note:
439: 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.
441: .seealso: [](chapter_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
442: @*/
443: PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
444: {
445: DMTS tdm;
446: DMTS_Local *dmlocalts;
447: const char *prefix;
450: DMGetDMTSWrite(dm, &tdm);
451: DMLocalTSGetContext(dm, tdm, &dmlocalts);
452: DMCreateMassMatrix(dm, dm, &dmlocalts->mass);
453: KSPCreate(PetscObjectComm((PetscObject)dm), &dmlocalts->kspmass);
454: PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix);
455: KSPSetOptionsPrefix(dmlocalts->kspmass, prefix);
456: KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_");
457: KSPSetFromOptions(dmlocalts->kspmass);
458: KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass);
459: return 0;
460: }
462: /*@C
463: 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.
465: Collective
467: Input Parameters:
468: . dm - `DM` providing the mass matrix
470: Level: developer
472: Note:
473: 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.
474: Since the matrix is lumped, inversion is trivial.
476: .seealso: [](chapter_ts), `DM`, `DMTSCreateRHSMassMatrix()`, `DMTSDestroyRHSMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
477: @*/
478: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
479: {
480: DMTS tdm;
481: DMTS_Local *dmlocalts;
484: DMGetDMTSWrite(dm, &tdm);
485: DMLocalTSGetContext(dm, tdm, &dmlocalts);
486: DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv);
487: VecReciprocal(dmlocalts->lumpedmassinv);
488: VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view");
489: return 0;
490: }
492: /*@C
493: DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the `DM` context, if they exist.
495: Logically Collective
497: Input Parameters:
498: . dm - `DM` providing the mass matrix
500: Level: developer
502: .seealso: [](chapter_ts), `DM`, `DMTSCreateRHSMassMatrixLumped()`, `DMCreateMassMatrix()`, `DMCreateMassMatrix()`, `DMTS`
503: @*/
504: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
505: {
506: DMTS tdm;
507: DMTS_Local *dmlocalts;
510: DMGetDMTSWrite(dm, &tdm);
511: DMLocalTSGetContext(dm, tdm, &dmlocalts);
512: VecDestroy(&dmlocalts->lumpedmassinv);
513: MatDestroy(&dmlocalts->mass);
514: KSPDestroy(&dmlocalts->kspmass);
515: return 0;
516: }