Actual source code: dmplexts.c
petsc-3.9.4 2018-09-11
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: PetscInt cStart, cEnd, cEndInterior;
109: DM plex;
113: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
114: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
115: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
116: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
117: DMGetLocalVector(plex, &locF);
118: VecZeroEntries(locF);
119: DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, NULL, time, locF, user);
120: DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
121: DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
122: DMRestoreLocalVector(plex, &locF);
123: DMDestroy(&plex);
124: return(0);
125: }
127: /*@
128: 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
130: Input Parameters:
131: + dm - The mesh
132: . t - The time
133: . locX - Local solution
134: . locX_t - Local solution time derivative, or NULL
135: - user - The user context
137: Level: developer
139: .seealso: DMPlexComputeJacobianActionFEM()
140: @*/
141: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
142: {
143: DM plex;
144: Vec faceGeometryFVM = NULL;
145: PetscInt Nf, f;
149: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
150: DMGetNumFields(plex, &Nf);
151: if (!locX_t) {
152: /* This is the RHS part */
153: for (f = 0; f < Nf; f++) {
154: PetscObject obj;
155: PetscClassId id;
157: DMGetField(plex, f, &obj);
158: PetscObjectGetClassId(obj, &id);
159: if (id == PETSCFV_CLASSID) {
160: DMPlexSNESGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
161: break;
162: }
163: }
164: }
165: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
166: /* TODO: locX_t */
167: DMDestroy(&plex);
168: return(0);
169: }
171: /*@
172: DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
174: Input Parameters:
175: + dm - The mesh
176: . t - The time
177: . locX - Local solution
178: . locX_t - Local solution time derivative, or NULL
179: - user - The user context
181: Output Parameter:
182: . locF - Local output vector
184: Level: developer
186: .seealso: DMPlexComputeJacobianActionFEM()
187: @*/
188: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
189: {
190: PetscInt cStart, cEnd, cEndInterior;
191: DM plex;
195: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
196: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
197: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
198: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
199: DMPlexComputeResidual_Internal(plex, cStart, cEnd, time, locX, locX_t, time, locF, user);
200: DMDestroy(&plex);
201: return(0);
202: }
204: /*@
205: DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
207: Input Parameters:
208: + dm - The mesh
209: . t - The time
210: . locX - Local solution
211: . locX_t - Local solution time derivative, or NULL
212: . X_tshift - The multiplicative parameter for dF/du_t
213: - user - The user context
215: Output Parameter:
216: . locF - Local output vector
218: Level: developer
220: .seealso: DMPlexComputeJacobianActionFEM()
221: @*/
222: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
223: {
224: PetscInt cStart, cEnd, cEndInterior;
225: DM plex;
229: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
230: DMPlexGetHeightStratum(plex, 0, &cStart, &cEnd);
231: DMPlexGetHybridBounds(plex, &cEndInterior, NULL, NULL, NULL);
232: cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
233: DMPlexComputeJacobian_Internal(plex, cStart, cEnd, time, X_tShift, locX, locX_t, Jac, JacP, user);
234: DMDestroy(&plex);
235: return(0);
236: }
238: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
239: {
240: DM dm;
241: SNES snes;
242: Vec sol;
243: PetscBool check;
247: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
248: if (!check) return(0);
249: VecDuplicate(u, &sol);
250: TSSetSolution(ts, sol);
251: TSGetDM(ts, &dm);
252: TSSetUp(ts);
253: TSGetSNES(ts, &snes);
254: SNESSetSolution(snes, sol);
255: DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
256: VecDestroy(&sol);
257: return(0);
258: }