Actual source code: dmlocalts.c
petsc-3.7.3 2016-08-01
1: #include <petsc/private/dmimpl.h>
2: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
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;
17: static PetscErrorCode DMTSDestroy_DMLocal(DMTS tdm)
18: {
22: PetscFree(tdm->data);
23: return(0);
24: }
28: static PetscErrorCode DMTSDuplicate_DMLocal(DMTS oldtdm, DMTS tdm)
29: {
33: PetscNewLog(tdm, (DMTS_Local **) &tdm->data);
34: if (oldtdm->data) {PetscMemcpy(tdm->data, oldtdm->data, sizeof(DMTS_Local));}
35: return(0);
36: }
40: static PetscErrorCode DMLocalTSGetContext(DM dm, DMTS tdm, DMTS_Local **dmlocalts)
41: {
45: *dmlocalts = NULL;
46: if (!tdm->data) {
47: PetscNewLog(dm, (DMTS_Local **) &tdm->data);
49: tdm->ops->destroy = DMTSDestroy_DMLocal;
50: tdm->ops->duplicate = DMTSDuplicate_DMLocal;
51: }
52: *dmlocalts = (DMTS_Local *) tdm->data;
53: return(0);
54: }
58: static PetscErrorCode TSComputeIFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, Vec F, void *ctx)
59: {
60: DM dm;
61: Vec locX, locX_t, locF;
62: DMTS_Local *dmlocalts = (DMTS_Local *) ctx;
70: TSGetDM(ts, &dm);
71: DMGetLocalVector(dm, &locX);
72: DMGetLocalVector(dm, &locX_t);
73: DMGetLocalVector(dm, &locF);
74: VecZeroEntries(locX);
75: VecZeroEntries(locX_t);
76: if (dmlocalts->boundarylocal) {(*dmlocalts->boundarylocal)(dm, time, locX, locX_t,dmlocalts->boundarylocalctx);}
77: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
78: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
79: DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
80: DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
81: VecZeroEntries(locF);
82: CHKMEMQ;
83: (*dmlocalts->ifunctionlocal)(dm, time, locX, locX_t, locF, dmlocalts->ifunctionlocalctx);
84: CHKMEMQ;
85: VecZeroEntries(F);
86: DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F);
87: DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F);
88: DMRestoreLocalVector(dm, &locX);
89: DMRestoreLocalVector(dm, &locX_t);
90: DMRestoreLocalVector(dm, &locF);
91: return(0);
92: }
96: static PetscErrorCode TSComputeRHSFunction_DMLocal(TS ts, PetscReal time, Vec X, Vec F, void *ctx)
97: {
98: DM dm;
99: Vec locX;
100: DMTS_Local *dmlocalts = (DMTS_Local *) ctx;
107: TSGetDM(ts, &dm);
108: DMGetLocalVector(dm, &locX);
109: VecZeroEntries(locX);
110: if (dmlocalts->boundarylocal) {(*dmlocalts->boundarylocal)(dm,time,locX,NULL,dmlocalts->boundarylocalctx);}
111: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
112: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
113: VecZeroEntries(F);
114: CHKMEMQ;
115: (*dmlocalts->rhsfunctionlocal)(dm, time, locX, F, dmlocalts->rhsfunctionlocalctx);
116: CHKMEMQ;
117: DMRestoreLocalVector(dm, &locX);
118: return(0);
119: }
123: static PetscErrorCode TSComputeIJacobian_DMLocal(TS ts, PetscReal time, Vec X, Vec X_t, PetscReal a, Mat A, Mat B, void *ctx)
124: {
125: DM dm;
126: Vec locX, locX_t;
127: DMTS_Local *dmlocalts = (DMTS_Local *) ctx;
131: TSGetDM(ts, &dm);
132: if (dmlocalts->ijacobianlocal) {
133: DMGetLocalVector(dm, &locX);
134: DMGetLocalVector(dm, &locX_t);
135: VecZeroEntries(locX);
136: if (dmlocalts->boundarylocal) {(*dmlocalts->boundarylocal)(dm,time,locX,locX_t,dmlocalts->boundarylocalctx);}
137: VecZeroEntries(locX_t);
138: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX);
139: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX);
140: DMGlobalToLocalBegin(dm, X_t, INSERT_VALUES, locX_t);
141: DMGlobalToLocalEnd(dm, X_t, INSERT_VALUES, locX_t);
142: CHKMEMQ;
143: (*dmlocalts->ijacobianlocal)(dm, time, locX, locX_t, a, A, B, dmlocalts->ijacobianlocalctx);
144: CHKMEMQ;
145: DMRestoreLocalVector(dm, &locX);
146: DMRestoreLocalVector(dm, &locX_t);
147: } else {
148: MatFDColoring fdcoloring;
149: PetscObjectQuery((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject *) &fdcoloring);
150: if (!fdcoloring) {
151: ISColoring coloring;
153: DMCreateColoring(dm, dm->coloringtype, &coloring);
154: MatFDColoringCreate(B, coloring, &fdcoloring);
155: ISColoringDestroy(&coloring);
156: switch (dm->coloringtype) {
157: case IS_COLORING_GLOBAL:
158: MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void)) TSComputeIFunction_DMLocal, dmlocalts);
159: break;
160: default: SETERRQ1(PetscObjectComm((PetscObject) ts), PETSC_ERR_SUP, "No support for coloring type '%s'", ISColoringTypes[dm->coloringtype]);
161: }
162: PetscObjectSetOptionsPrefix((PetscObject) fdcoloring, ((PetscObject) dm)->prefix);
163: MatFDColoringSetFromOptions(fdcoloring);
164: MatFDColoringSetUp(B, coloring, fdcoloring);
165: PetscObjectCompose((PetscObject) dm, "DMDASNES_FDCOLORING", (PetscObject) fdcoloring);
166: PetscObjectDereference((PetscObject) fdcoloring);
168: /* The following breaks an ugly reference counting loop that deserves a paragraph. MatFDColoringApply() will call
169: * VecDuplicate() with the state Vec and store inside the MatFDColoring. This Vec will duplicate the Vec, but the
170: * MatFDColoring is composed with the DM. We dereference the DM here so that the reference count will eventually
171: * drop to 0. Note the code in DMDestroy() that exits early for a negative reference count. That code path will be
172: * taken when the PetscObjectList for the Vec inside MatFDColoring is destroyed.
173: */
174: PetscObjectDereference((PetscObject) dm);
175: }
176: MatFDColoringApply(B, fdcoloring, X, ts);
177: }
178: /* This will be redundant if the user called both, but it's too common to forget. */
179: if (A != B) {
180: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
181: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
182: }
183: return(0);
184: }
188: /*@C
189: DMTSSetBoundaryLocal - set the function for essential boundary data for a local implicit function evaluation.
190: 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).
191: Vectors are initialized to zero before this function, so it is only needed for non homogeneous data.
193: Note that this function is somewhat optional: boundary data could potentially be inserted by a function passed to
194: DMTSSetIFunctionLocal(). The use case for this function is for discretizations with constraints (see
195: DMGetDefaultConstraints()): this function inserts boundary values before constraint interpolation.
197: Logically Collective
199: Input Arguments:
200: + dm - DM to associate callback with
201: . func - local function evaluation
202: - ctx - context for function evaluation
204: Level: intermediate
206: .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal()
207: @*/
208: PetscErrorCode DMTSSetBoundaryLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
209: {
210: DMTS tdm;
211: DMTS_Local *dmlocalts;
216: DMGetDMTSWrite(dm, &tdm);
217: DMLocalTSGetContext(dm, tdm, &dmlocalts);
219: dmlocalts->boundarylocal = func;
220: dmlocalts->boundarylocalctx = ctx;
222: return(0);
223: }
227: /*@C
228: DMTSSetIFunctionLocal - set a local implicit function evaluation function. This function is called with local vector
229: containing the local vector information PLUS ghost point information. It should compute a result for all local
230: elements and DMTS will automatically accumulate the overlapping values.
232: Logically Collective
234: Input Arguments:
235: + dm - DM to associate callback with
236: . func - local function evaluation
237: - ctx - context for function evaluation
239: Level: beginner
241: .seealso: DMTSSetIFunction(), DMTSSetIJacobianLocal()
242: @*/
243: PetscErrorCode DMTSSetIFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, Vec, void *), void *ctx)
244: {
245: DMTS tdm;
246: DMTS_Local *dmlocalts;
251: DMGetDMTSWrite(dm, &tdm);
252: DMLocalTSGetContext(dm, tdm, &dmlocalts);
254: dmlocalts->ifunctionlocal = func;
255: dmlocalts->ifunctionlocalctx = ctx;
257: DMTSSetIFunction(dm, TSComputeIFunction_DMLocal, dmlocalts);
258: if (!tdm->ops->ijacobian) { /* Call us for the Jacobian too, can be overridden by the user. */
259: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
260: }
261: return(0);
262: }
266: /*@C
267: DMTSSetIJacobianLocal - set a local Jacobian evaluation function
269: Logically Collective
271: Input Arguments:
272: + dm - DM to associate callback with
273: . func - local Jacobian evaluation
274: - ctx - optional context for local Jacobian evaluation
276: Level: beginner
278: .seealso: DMTSSetIFunctionLocal(), DMTSSetIJacobian(), DMTSSetIFunction()
279: @*/
280: PetscErrorCode DMTSSetIJacobianLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *), void *ctx)
281: {
282: DMTS tdm;
283: DMTS_Local *dmlocalts;
288: DMGetDMTSWrite(dm, &tdm);
289: DMLocalTSGetContext(dm, tdm, &dmlocalts);
291: dmlocalts->ijacobianlocal = func;
292: dmlocalts->ijacobianlocalctx = ctx;
294: DMTSSetIJacobian(dm, TSComputeIJacobian_DMLocal, dmlocalts);
295: return(0);
296: }
300: /*@C
301: DMTSSetRHSFunctionLocal - set a local rhs function evaluation function. This function is called with local vector
302: containing the local vector information PLUS ghost point information. It should compute a result for all local
303: elements and DMTS will automatically accumulate the overlapping values.
305: Logically Collective
307: Input Arguments:
308: + dm - DM to associate callback with
309: . func - local function evaluation
310: - ctx - context for function evaluation
312: Level: beginner
314: .seealso: DMTSSetRHSFunction(), DMTSSetIFunction(), DMTSSetIJacobianLocal()
315: @*/
316: PetscErrorCode DMTSSetRHSFunctionLocal(DM dm, PetscErrorCode (*func)(DM, PetscReal, Vec, Vec, void *), void *ctx)
317: {
318: DMTS tdm;
319: DMTS_Local *dmlocalts;
324: DMGetDMTSWrite(dm, &tdm);
325: DMLocalTSGetContext(dm, tdm, &dmlocalts);
327: dmlocalts->rhsfunctionlocal = func;
328: dmlocalts->rhsfunctionlocalctx = ctx;
330: DMTSSetRHSFunction(dm, TSComputeRHSFunction_DMLocal, dmlocalts);
331: return(0);
332: }