Actual source code: dmlocalts.c
petsc-3.14.6 2021-03-30
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: } DMTS_Local;
15: static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm)
16: {
20: PetscFree(tdm->data);
21: return(0);
22: }
24: static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
25: {
29: PetscNewLog(tdm, (DMTS_Local **) &tdm->data);
30: if (oldtdm->data) {PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local));}
31: return(0);
32: }
34: static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts)
35: {
39: *dmlocalts = NULL;
40: if (!tdm->data) {
41: PetscNewLog(dm, (DMTS_Local **) &tdm->data);
43: tdm->ops->destroy = DMTSDestroy_DMLocal;
44: tdm->ops->duplicate = DMTSDuplicate_DMLocal;
45: }
46: *dmlocalts = (DMTS_Local *) tdm->data;
47: return(0);
48: }
50: static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, void *ctx)
51: {
52: DM dm;
53: Vec locX, locX_t, locF;
54: DMTS_Local *dmlocalts = (DMTS_Local *) ctx;
62: TSGetDM(ts, &dm);
63: DMGetLocalVector(dm, &locX);
64: DMGetLocalVector(dm, &locX_t);
65: DMGetLocalVector(dm, &locF);
66: VecZeroEntries(locX);
67: VecZeroEntries(locX_t);
68: if (dmlocalts->boundarylocal) {(*dmlocalts->boundarylocal)(dm, time, locX, locX_t,dmlocalts->boundarylocalctx);}
69: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
70: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
71: DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
72: DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
73: VecZeroEntries(locF);
74: CHKMEMQ;
75: (*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx);
76: CHKMEMQ;
77: VecZeroEntries(F);
78: DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F);
79: DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F);
80: DMRestoreLocalVector(dm, &locX);
81: DMRestoreLocalVector(dm, &locX_t);
82: DMRestoreLocalVector(dm, &locF);
83: return(0);
84: }
86: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, void *ctx)
87: {
88: DM dm;
89: Vec locX;
90: DMTS_Local *dmlocalts = (DMTS_Local *) ctx;
97: TSGetDM(ts, &dm);
98: DMGetLocalVector(dm, &locX);
99: VecZeroEntries(locX);
100: if (dmlocalts->boundarylocal) {(*dmlocalts->boundarylocal)(dm,time,locX,NULL,dmlocalts->boundarylocalctx);}
101: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
102: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
103: VecZeroEntries(F);
104: CHKMEMQ;
105: (*dmlocalts->rhsfunctionlocal)(dm, time, locX, F, dmlocalts->rhsfunctionlocalctx);
106: CHKMEMQ;
107: DMRestoreLocalVector(dm, &locX);
108: return(0);
109: }
111: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, void *ctx)
112: {
113: DM dm;
114: Vec locX, locX_t;
115: DMTS_Local *dmlocalts = (DMTS_Local *) ctx;
119: TSGetDM(ts, &dm);
120: if (dmlocalts->ijacobianlocal) {
121: DMGetLocalVector(dm, &locX);
122: DMGetLocalVector(dm, &locX_t);
123: VecZeroEntries(locX);
124: VecZeroEntries(locX_t);
125: if (dmlocalts->boundarylocal) {(*dmlocalts->boundarylocal)(dm,time,locX,locX_t,dmlocalts->boundarylocalctx);}
126: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
127: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
128: DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
129: DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
130: CHKMEMQ;
131: (*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx);
132: CHKMEMQ;
133: DMRestoreLocalVector(dm, &locX);
134: DMRestoreLocalVector(dm, &locX_t);
135: } else {
136: MatFDColoring fdcoloring;
137: PetscObjectQuery((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject *) &fdcoloring);
138: if (!fdcoloring) {
139: ISColoring coloring;
141: DMCreateColoring(dm, dm->coloringtype, &coloring);
142: MatFDColoringCreate(B, coloring, &fdcoloring);
143: ISColoringDestroy(&coloring);
144: switch (dm->coloringtype) {
145: case IS_COLORING_GLOBAL:
146: MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void)) TSComputeIFunction_DMLocal, dmlocalts);
147: break;
148: default: SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
149: }
150: PetscObjectSetOptionsPrefix((PetscObject) fdcoloring, ((PetscObject) dm)->prefix);
151: MatFDColoringSetFromOptions(fdcoloring);
152: MatFDColoringSetUp(B, coloring, fdcoloring);
153: PetscObjectCompose((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject) fdcoloring);
154: PetscObjectDereference((PetscObject) fdcoloring);
156: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
157: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
158: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
159: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
160: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
161: */
162: PetscObjectDereference((PetscObject) dm);
163: }
164: MatFDColoringApply(B, fdcoloring, X, ts);
165: }
166: /* This will be redundant if the user called both, but it's too common to forget. */
167: if (A != B) {
168: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
169: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
170: }
171: return(0);
172: }
174: /*@C
175: DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
176: 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).
177: Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.
179: Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to
180: DMTSSetIFunctionLocal(). The use case for this function is for discretizations with constraints (see
181: DMGetDefaultConstraints()): this function inserts boundary values before constraint interpolation.
183: Logically Collective
185: Input Arguments:
186: + dm - DM to associate callback with
187: . func - local function evaluation
188: - ctx - context for function evaluation
190: Level: intermediate
192: .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal()
193: @*/
194: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
195: {
196: DMTS tdm;
197: DMTS_Local *dmlocalts;
202: DMGetDMTSWrite(dm, &tdm);
203: DMLocalTSGetContext(dm, tdm, &dmlocalts);
205: dmlocalts->boundarylocal = func;
206: dmlocalts->boundarylocalctx = ctx;
208: return(0);
209: }
211: /*@C
212: DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
213: containing the local vector information PLUS ghost point information. It should compute a result for all local
214: elements and DMTS will automatically accumulate the overlapping values.
216: Logically Collective
218: Input Arguments:
219: + dm - DM to associate callback with
220: . func - local function evaluation
221: - ctx - context for function evaluation
223: Level: beginner
225: .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal()
226: @*/
227: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
228: {
229: DMTS tdm;
230: DMTS_Local *dmlocalts;
235: DMGetDMTSWrite(dm, &tdm);
236: DMLocalTSGetContext(dm, tdm, &dmlocalts);
238: dmlocalts->ifunctionlocal = func;
239: dmlocalts->ifunctionlocalctx = ctx;
241: DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts);
242: if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
243: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
244: }
245: return(0);
246: }
248: /*@C
249: DMTSSetIJacobianLocal - set a local Jacobian evaluation function
251: Logically Collective
253: Input Arguments:
254: + dm - DM to associate callback with
255: . func - local Jacobian evaluation
256: - ctx - optional context for local Jacobian evaluation
258: Level: beginner
260: .seealso: DMTSSetIFunctionLocal(), DMTSSetIJacobian(), DMTSSetIFunction()
261: @*/
262: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
263: {
264: DMTS tdm;
265: DMTS_Local *dmlocalts;
270: DMGetDMTSWrite(dm, &tdm);
271: DMLocalTSGetContext(dm, tdm, &dmlocalts);
273: dmlocalts->ijacobianlocal = func;
274: dmlocalts->ijacobianlocalctx = ctx;
276: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
277: return(0);
278: }
280: /*@C
281: DMTSSetRHSFunctionLocal - set a local rhs 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 DMTS will automatically accumulate the overlapping values.
285: Logically Collective
287: Input Arguments:
288: + dm - DM to associate callback with
289: . func - local function evaluation
290: - ctx - context for function evaluation
292: Level: beginner
294: .seealso: DMTSSetRHSFunction(), DMTSSetIFunction(), DMTSSetIJacobianLocal()
295: @*/
296: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
297: {
298: DMTS tdm;
299: DMTS_Local *dmlocalts;
304: DMGetDMTSWrite(dm, &tdm);
305: DMLocalTSGetContext(dm, tdm, &dmlocalts);
307: dmlocalts->rhsfunctionlocal = func;
308: dmlocalts->rhsfunctionlocalctx = ctx;
310: DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts);
311: return(0);
312: }