Actual source code: dmplexts.c
petsc-3.11.4 2019-09-28
1: #include <petsc/private/dmpleximpl.h>
2: #include <petsc/private/tsimpl.h>
3: #include <petsc/private/snesimpl.h>
4: #include <petscds.h>
5: #include <petscfv.h>
7: static PetscErrorCode DMTSConvertPlex(DM dm, DM *plex, PetscBool copy)
8: {
9: PetscBool isPlex;
13: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
14: if (isPlex) {
15: *plex = dm;
16: PetscObjectReference((PetscObject) dm);
17: } else {
18: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
19: if (!*plex) {
20: DMConvert(dm,DMPLEX,plex);
21: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
22: if (copy) {
23: PetscInt i;
24: PetscObject obj;
25: const char *comps[3] = {"A","dmAux","dmCh"};
27: DMCopyDMTS(dm, *plex);
28: DMCopyDMSNES(dm, *plex);
29: for (i = 0; i < 3; i++) {
30: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
31: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
32: }
33: }
34: } else {
35: PetscObjectReference((PetscObject) *plex);
36: }
37: }
38: return(0);
39: }
42: /*@
43: DMPlexTSGetGeometryFVM - Return precomputed geometric data
45: Input Parameter:
46: . dm - The DM
48: Output Parameters:
49: + facegeom - The values precomputed from face geometry
50: . cellgeom - The values precomputed from cell geometry
51: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell
53: Level: developer
55: .seealso: DMPlexTSSetRHSFunctionLocal()
56: @*/
57: PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
58: {
62: DMPlexSNESGetGeometryFVM(dm,facegeom,cellgeom,minRadius);
63: return(0);
64: }
66: /*@C
67: DMPlexTSGetGradientDM - Return gradient data layout
69: Input Parameters:
70: + dm - The DM
71: - fv - The PetscFV
73: Output Parameter:
74: . dmGrad - The layout for gradient values
76: Level: developer
78: .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal()
79: @*/
80: PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
81: {
85: DMPlexSNESGetGradientDM(dm,fv,dmGrad);
86: return(0);
87: }
89: /*@
90: DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user
92: Input Parameters:
93: + dm - The mesh
94: . t - The time
95: . locX - Local solution
96: - user - The user context
98: Output Parameter:
99: . F - Global output vector
101: Level: developer
103: .seealso: DMPlexComputeJacobianActionFEM()
104: @*/
105: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
106: {
107: Vec locF;
108: IS cellIS;
109: DM plex;
110: PetscInt depth;
114: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
115: DMPlexGetDepth(plex, &depth);
116: DMGetStratumIS(plex, "dim", depth, &cellIS);
117: if (!cellIS) {
118: DMGetStratumIS(plex, "depth", depth, &cellIS);
119: }
120: DMGetLocalVector(plex, &locF);
121: VecZeroEntries(locF);
122: DMPlexComputeResidual_Internal(plex, cellIS, time, locX, NULL, time, locF, user);
123: DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
124: DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
125: DMRestoreLocalVector(plex, &locF);
126: ISDestroy(&cellIS);
127: DMDestroy(&plex);
128: return(0);
129: }
131: /*@
132: 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
134: Input Parameters:
135: + dm - The mesh
136: . t - The time
137: . locX - Local solution
138: . locX_t - Local solution time derivative, or NULL
139: - user - The user context
141: Level: developer
143: .seealso: DMPlexComputeJacobianActionFEM()
144: @*/
145: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
146: {
147: DM plex;
148: Vec faceGeometryFVM = NULL;
149: PetscInt Nf, f;
153: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
154: DMGetNumFields(plex, &Nf);
155: if (!locX_t) {
156: /* This is the RHS part */
157: for (f = 0; f < Nf; f++) {
158: PetscObject obj;
159: PetscClassId id;
161: DMGetField(plex, f, NULL, &obj);
162: PetscObjectGetClassId(obj, &id);
163: if (id == PETSCFV_CLASSID) {
164: DMPlexSNESGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
165: break;
166: }
167: }
168: }
169: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
170: /* TODO: locX_t */
171: DMDestroy(&plex);
172: return(0);
173: }
175: /*@
176: DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
178: Input Parameters:
179: + dm - The mesh
180: . t - The time
181: . locX - Local solution
182: . locX_t - Local solution time derivative, or NULL
183: - user - The user context
185: Output Parameter:
186: . locF - Local output vector
188: Level: developer
190: .seealso: DMPlexComputeJacobianActionFEM()
191: @*/
192: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
193: {
194: DM plex;
195: IS cellIS;
196: PetscInt depth;
200: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
201: DMPlexGetDepth(plex, &depth);
202: DMGetStratumIS(plex, "dim", depth, &cellIS);
203: if (!cellIS) {
204: DMGetStratumIS(plex, "depth", depth, &cellIS);
205: }
206: DMPlexComputeResidual_Internal(plex, cellIS, time, locX, locX_t, time, locF, user);
207: ISDestroy(&cellIS);
208: DMDestroy(&plex);
209: return(0);
210: }
212: /*@
213: DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
215: Input Parameters:
216: + dm - The mesh
217: . t - The time
218: . locX - Local solution
219: . locX_t - Local solution time derivative, or NULL
220: . X_tshift - The multiplicative parameter for dF/du_t
221: - user - The user context
223: Output Parameter:
224: . locF - Local output vector
226: Level: developer
228: .seealso: DMPlexComputeJacobianActionFEM()
229: @*/
230: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
231: {
232: DM plex;
233: PetscDS prob;
234: PetscBool hasJac, hasPrec;
235: IS cellIS;
236: PetscInt depth;
240: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
241: DMGetDS(dm, &prob);
242: PetscDSHasJacobian(prob, &hasJac);
243: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
244: if (hasJac && hasPrec) {MatZeroEntries(Jac);}
245: MatZeroEntries(JacP);
246: DMPlexGetDepth(plex,&depth);
247: DMGetStratumIS(plex, "dim", depth, &cellIS);
248: if (!cellIS) {DMGetStratumIS(plex, "depth", depth, &cellIS);}
249: DMPlexComputeJacobian_Internal(plex, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
250: ISDestroy(&cellIS);
251: DMDestroy(&plex);
252: return(0);
253: }
255: /*@C
256: DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
258: Input Parameters:
259: + ts - the TS object
260: . u - representative TS vector
261: . exactFuncs - pointwise functions of the exact solution for each field
262: - ctxs - contexts for the functions
264: Level: developer
265: @*/
266: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
267: {
268: DM dm;
269: SNES snes;
270: Vec sol;
271: PetscBool check;
275: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
276: if (!check) return(0);
277: VecDuplicate(u, &sol);
278: VecCopy(u, sol);
279: TSSetSolution(ts, u);
280: TSGetDM(ts, &dm);
281: TSSetUp(ts);
282: TSGetSNES(ts, &snes);
283: SNESSetSolution(snes, u);
284: DMSNESCheckFromOptions_Internal(snes, dm, sol, exactFuncs, ctxs);
285: VecDestroy(&sol);
286: return(0);
287: }