Actual source code: dmplexts.c
petsc-3.7.5 2017-01-01
1: #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
2: #include <petsc/private/tsimpl.h> /*I "petscts.h" I*/
3: #include <petsc/private/snesimpl.h>
4: #include <petscds.h>
5: #include <petscfv.h>
9: static PetscErrorCode DMTSConvertPlex(DM dm, DM *plex, PetscBool copy)
10: {
11: PetscBool isPlex;
15: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
16: if (isPlex) {
17: *plex = dm;
18: PetscObjectReference((PetscObject) dm);
19: } else {
20: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
21: if (!*plex) {
22: DMConvert(dm,DMPLEX,plex);
23: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
24: if (copy) {
25: PetscInt i;
26: PetscObject obj;
27: const char *comps[3] = {"A","dmAux","dmCh"};
29: DMCopyDMTS(dm, *plex);
30: DMCopyDMSNES(dm, *plex);
31: for (i = 0; i < 3; i++) {
32: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
33: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
34: }
35: }
36: } else {
37: PetscObjectReference((PetscObject) *plex);
38: }
39: }
40: return(0);
41: }
46: /*@
47: DMPlexTSGetGeometryFVM - Return precomputed geometric data
49: Input Parameter:
50: . dm - The DM
52: Output Parameters:
53: + facegeom - The values precomputed from face geometry
54: . cellgeom - The values precomputed from cell geometry
55: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
57: Level: developer
59: .seealso: DMPlexTSSetRHSFunctionLocal()
60: @*/
61: PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
62: {
63: DMTS dmts;
64: PetscObject obj;
65: DM plex;
70: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
71: DMGetDMTS(plex, &dmts);
72: PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", &obj);
73: if (!obj) {
74: Vec cellgeom, facegeom;
76: DMPlexComputeGeometryFVM(plex, &cellgeom, &facegeom);
77: PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject) facegeom);
78: PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject) cellgeom);
79: VecDestroy(&facegeom);
80: VecDestroy(&cellgeom);
81: }
84: if (minRadius) {DMPlexGetMinRadius(plex, minRadius);}
85: DMDestroy(&plex);
86: return(0);
87: }
91: /*@C
92: DMPlexTSGetGradientDM - Return gradient data layout
94: Input Parameters:
95: + dm - The DM
96: - fv - The PetscFV
98: Output Parameter:
99: . dmGrad - The layout for gradient values
101: Level: developer
103: .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal()
104: @*/
105: PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
106: {
107: DMTS dmts;
108: PetscObject obj;
109: PetscBool computeGradients;
110: DM plex;
117: PetscFVGetComputeGradients(fv, &computeGradients);
118: if (!computeGradients) {*dmGrad = NULL; return(0);}
119: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
120: DMGetDMTS(plex, &dmts);
121: PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);
122: if (!obj) {
123: DM dmGrad;
124: Vec faceGeometry, cellGeometry;
126: DMPlexTSGetGeometryFVM(plex, &faceGeometry, &cellGeometry, NULL);
127: DMPlexComputeGradientFVM(plex, fv, faceGeometry, cellGeometry, &dmGrad);
128: PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);
129: DMDestroy(&dmGrad);
130: }
131: PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);
132: DMDestroy(&plex);
133: return(0);
134: }
138: /*@
139: DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user
141: Input Parameters:
142: + dm - The mesh
143: . t - The time
144: . locX - Local solution
145: - user - The user context
147: Output Parameter:
148: . F - Global output vector
150: Level: developer
152: .seealso: DMPlexComputeJacobianActionFEM()
153: @*/
154: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
155: {
156: Vec locF;
157: PetscInt cStart, cEnd, cEndInterior;
158: DM plex;
162: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
163: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
164: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
165: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
166: DMGetLocalVector(plex, &locF);
167: VecZeroEntries(locF);
168: DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, NULL, locF, user);
169: DMLocalToGlobalBegin(plex, locF, INSERT_VALUES, F);
170: DMLocalToGlobalEnd(plex, locF, INSERT_VALUES, F);
171: DMRestoreLocalVector(plex, &locF);
172: DMDestroy(&plex);
173: return(0);
174: }
178: /*@
179: DMPlexTSComputeBoundary - Insert the essential boundary values for the local input X and/or its time derivative X_t using pointwise functions specified by the user
181: Input Parameters:
182: + dm - The mesh
183: . t - The time
184: . locX - Local solution
185: . locX_t - Local solution time derivative, or NULL
186: - user - The user context
188: Level: developer
190: .seealso: DMPlexComputeJacobianActionFEM()
191: @*/
192: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
193: {
194: DM plex;
195: Vec faceGeometryFVM = NULL;
196: PetscInt Nf, f;
200: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
201: DMGetNumFields(plex, &Nf);
202: if (!locX_t) {
203: /* This is the RHS part */
204: for (f = 0; f < Nf; f++) {
205: PetscObject obj;
206: PetscClassId id;
208: DMGetField(plex, f, &obj);
209: PetscObjectGetClassId(obj, &id);
210: if (id == PETSCFV_CLASSID) {
211: DMPlexSNESGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
212: break;
213: }
214: }
215: }
216: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
217: /* TODO: locX_t */
218: DMDestroy(&plex);
219: return(0);
220: }
224: /*@
225: DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
227: Input Parameters:
228: + dm - The mesh
229: . t - The time
230: . locX - Local solution
231: . locX_t - Local solution time derivative, or NULL
232: - user - The user context
234: Output Parameter:
235: . locF - Local output vector
237: Level: developer
239: .seealso: DMPlexComputeJacobianActionFEM()
240: @*/
241: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
242: {
243: PetscInt cStart, cEnd, cEndInterior;
244: DM plex;
248: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
249: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
250: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
251: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
252: DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, locX_t, locF, user);
253: DMDestroy(&plex);
254: return(0);
255: }
259: /*@
260: DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
262: Input Parameters:
263: + dm - The mesh
264: . t - The time
265: . locX - Local solution
266: . locX_t - Local solution time derivative, or NULL
267: . X_tshift - The multiplicative parameter for dF/du_t
268: - user - The user context
270: Output Parameter:
271: . locF - Local output vector
273: Level: developer
275: .seealso: DMPlexComputeJacobianActionFEM()
276: @*/
277: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
278: {
279: PetscInt cStart, cEnd, cEndInterior;
280: DM plex;
284: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
285: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
286: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
287: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
288: DMPlexComputeJacobian_Internal(plex, cStart, cEnd, time, X_tShift, locX, locX_t, Jac, JacP, user);
289: DMDestroy(&plex);
290: return(0);
291: }
295: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
296: {
297: DM dm;
298: SNES snes;
299: Vec sol;
300: PetscBool check;
304: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
305: if (!check) return(0);
306: VecDuplicate(u, &sol);
307: TSSetSolution(ts, sol);
308: TSGetDM(ts, &dm);
309: TSSetUp(ts);
310: TSGetSNES(ts, &snes);
311: SNESSetSolution(snes, sol);
312: DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
313: VecDestroy(&sol);
314: return(0);
315: }