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: PetscNewLog(tdm, (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: PetscNewLog(dm, (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: SETERRQ(PetscObjectComm((PetscObject) ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
152: }
153: PetscObjectSetOptionsPrefix((PetscObject) fdcoloring, ((PetscObject) dm)->prefix);
154: MatFDColoringSetFromOptions(fdcoloring);
155: MatFDColoringSetUp(B, coloring, fdcoloring);
156: PetscObjectCompose((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject) fdcoloring);
157: PetscObjectDereference((PetscObject) fdcoloring);
159: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
160: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
161: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
162: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
163: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
164: */
165: PetscObjectDereference((PetscObject) dm);
166: }
167: MatFDColoringApply(B, fdcoloring, X, ts);
168: }
169: /* This will be redundant if the user called both, but it's too common to forget. */
170: if (A != B) {
171: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
172: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
173: }
174: return 0;
175: }
177: /*@C
178: DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
179: 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).
180: Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.
182: Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to
183: DMTSSetIFunctionLocal(). The use case for this function is for discretizations with constraints (see
184: DMGetDefaultConstraints()): this function inserts boundary values before constraint interpolation.
186: Logically Collective
188: Input Parameters:
189: + dm - DM to associate callback with
190: . func - local function evaluation
191: - ctx - context for function evaluation
193: Level: intermediate
195: .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal()
196: @*/
197: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
198: {
199: DMTS tdm;
200: DMTS_Local *dmlocalts;
203: DMGetDMTSWrite(dm, &tdm);
204: DMLocalTSGetContext(dm, tdm, &dmlocalts);
206: dmlocalts->boundarylocal = func;
207: dmlocalts->boundarylocalctx = ctx;
209: return 0;
210: }
212: /*@C
213: DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
214: containing the local vector information PLUS ghost point information. It should compute a result for all local
215: elements and DMTS will automatically accumulate the overlapping values.
217: Logically Collective
219: Input Parameters:
220: + dm - DM to associate callback with
221: . func - local function evaluation
222: - ctx - context for function evaluation
224: Level: beginner
226: .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal()
227: @*/
228: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
229: {
230: DMTS tdm;
231: DMTS_Local *dmlocalts;
234: DMGetDMTSWrite(dm, &tdm);
235: DMLocalTSGetContext(dm, tdm, &dmlocalts);
237: dmlocalts->ifunctionlocal = func;
238: dmlocalts->ifunctionlocalctx = ctx;
240: DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts);
241: if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
242: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
243: }
244: return 0;
245: }
247: /*@C
248: DMTSSetIJacobianLocal - set a local Jacobian evaluation function
250: Logically Collective
252: Input Parameters:
253: + dm - DM to associate callback with
254: . func - local Jacobian evaluation
255: - ctx - optional context for local Jacobian evaluation
257: Level: beginner
259: .seealso: DMTSSetIFunctionLocal(), DMTSSetIJacobian(), DMTSSetIFunction()
260: @*/
261: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
262: {
263: DMTS tdm;
264: DMTS_Local *dmlocalts;
267: DMGetDMTSWrite(dm, &tdm);
268: DMLocalTSGetContext(dm, tdm, &dmlocalts);
270: dmlocalts->ijacobianlocal = func;
271: dmlocalts->ijacobianlocalctx = ctx;
273: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
274: return 0;
275: }
277: /*@C
278: DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
279: containing the local vector information PLUS ghost point information. It should compute a result for all local
280: elements and DMTS will automatically accumulate the overlapping values.
282: Logically Collective
284: Input Parameters:
285: + dm - DM to associate callback with
286: . func - local function evaluation
287: - ctx - context for function evaluation
289: Level: beginner
291: .seealso: DMTSSetRHSFunction(), DMTSSetIFunction(), DMTSSetIJacobianLocal()
292: @*/
293: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
294: {
295: DMTS tdm;
296: DMTS_Local *dmlocalts;
299: DMGetDMTSWrite(dm, &tdm);
300: DMLocalTSGetContext(dm, tdm, &dmlocalts);
302: dmlocalts->rhsfunctionlocal = func;
303: dmlocalts->rhsfunctionlocalctx = ctx;
305: DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts);
306: return 0;
307: }
309: /*@C
310: DMTSCreateRHSMassMatrix - This creates the mass matrix associated with the given DM, and a solver to invert it, and stores them in the DMTS context.
312: Collective on dm
314: Input Parameters:
315: . dm - DM providing the mass matrix
317: Note: 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.
319: Level: developer
321: .seealso: DMTSCreateRHSMassMatrixLumped(), DMTSDestroyRHSMassMatrix(), DMCreateMassMatrix(), DMTS
322: @*/
323: PetscErrorCode DMTSCreateRHSMassMatrix(DM dm)
324: {
325: DMTS tdm;
326: DMTS_Local *dmlocalts;
327: const char *prefix;
330: DMGetDMTSWrite(dm, &tdm);
331: DMLocalTSGetContext(dm, tdm, &dmlocalts);
332: DMCreateMassMatrix(dm, dm, &dmlocalts->mass);
333: KSPCreate(PetscObjectComm((PetscObject) dm), &dmlocalts->kspmass);
334: PetscObjectGetOptionsPrefix((PetscObject) dm, &prefix);
335: KSPSetOptionsPrefix(dmlocalts->kspmass, prefix);
336: KSPAppendOptionsPrefix(dmlocalts->kspmass, "mass_");
337: KSPSetFromOptions(dmlocalts->kspmass);
338: KSPSetOperators(dmlocalts->kspmass, dmlocalts->mass, dmlocalts->mass);
339: return 0;
340: }
342: /*@C
343: DMTSCreateRHSMassMatrixLumped - This creates the lumped mass matrix associated with the given DM, and a solver to invert it, and stores them in the DMTS context.
345: Collective on dm
347: Input Parameters:
348: . dm - DM providing the mass matrix
350: Note: 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.
351: Since the matrix is lumped, inversion is trivial.
353: Level: developer
355: .seealso: DMTSCreateRHSMassMatrix(), DMTSDestroyRHSMassMatrix(), DMCreateMassMatrix(), DMTS
356: @*/
357: PetscErrorCode DMTSCreateRHSMassMatrixLumped(DM dm)
358: {
359: DMTS tdm;
360: DMTS_Local *dmlocalts;
363: DMGetDMTSWrite(dm, &tdm);
364: DMLocalTSGetContext(dm, tdm, &dmlocalts);
365: DMCreateMassMatrixLumped(dm, &dmlocalts->lumpedmassinv);
366: VecReciprocal(dmlocalts->lumpedmassinv);
367: VecViewFromOptions(dmlocalts->lumpedmassinv, NULL, "-lumped_mass_inv_view");
368: return 0;
369: }
371: /*@C
372: DMTSDestroyRHSMassMatrix - Destroys the mass matrix and solver stored in the DMTS context, if they exist.
374: Logically Collective
376: Input Parameters:
377: . dm - DM providing the mass matrix
379: Level: developer
381: .seealso: DMTSCreateRHSMassMatrixLumped(), DMCreateMassMatrix(), DMCreateMassMatrix(), DMTS
382: @*/
383: PetscErrorCode DMTSDestroyRHSMassMatrix(DM dm)
384: {
385: DMTS tdm;
386: DMTS_Local *dmlocalts;
389: DMGetDMTSWrite(dm, &tdm);
390: DMLocalTSGetContext(dm, tdm, &dmlocalts);
391: VecDestroy(&dmlocalts->lumpedmassinv);
392: MatDestroy(&dmlocalts->mass);
393: KSPDestroy(&dmlocalts->kspmass);
394: return 0;
395: }