Actual source code: dmplexts.c

petsc-3.7.7 2017-09-25
Report Typos and Errors
  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: }