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: /*@
 10:   DMPlexTSGetGeometryFVM - Return precomputed geometric data

 12:   Input Parameter:
 13: . dm - The DM

 15:   Output Parameters:
 16: + facegeom - The values precomputed from face geometry
 17: . cellgeom - The values precomputed from cell geometry
 18: - minRadius - The minimum radius over the mesh of an inscribed sphere in a cell

 20:   Level: developer

 22: .seealso: DMPlexTSSetRHSFunctionLocal()
 23: @*/
 24: PetscErrorCode DMPlexTSGetGeometryFVM(DM dm, Vec *facegeom, Vec *cellgeom, PetscReal *minRadius)
 25: {
 26:   DMTS           dmts;
 27:   PetscObject    obj;

 32:   DMGetDMTS(dm, &dmts);
 33:   PetscObjectQuery((PetscObject) dmts, "DMPlexTS_facegeom_fvm", &obj);
 34:   if (!obj) {
 35:     Vec cellgeom, facegeom;

 37:     DMPlexComputeGeometryFVM(dm, &cellgeom, &facegeom);
 38:     PetscObjectCompose((PetscObject) dmts, "DMPlexTS_facegeom_fvm", (PetscObject) facegeom);
 39:     PetscObjectCompose((PetscObject) dmts, "DMPlexTS_cellgeom_fvm", (PetscObject) cellgeom);
 40:     VecDestroy(&facegeom);
 41:     VecDestroy(&cellgeom);
 42:   }
 45:   if (minRadius) {DMPlexGetMinRadius(dm, minRadius);}
 46:   return(0);
 47: }

 51: /*@C
 52:   DMPlexTSGetGradientDM - Return gradient data layout

 54:   Input Parameters:
 55: + dm - The DM
 56: - fv - The PetscFV

 58:   Output Parameter:
 59: . dmGrad - The layout for gradient values

 61:   Level: developer

 63: .seealso: DMPlexTSGetGeometryFVM(), DMPlexTSSetRHSFunctionLocal()
 64: @*/
 65: PetscErrorCode DMPlexTSGetGradientDM(DM dm, PetscFV fv, DM *dmGrad)
 66: {
 67:   DMTS           dmts;
 68:   PetscObject    obj;
 69:   PetscBool      computeGradients;

 76:   PetscFVGetComputeGradients(fv, &computeGradients);
 77:   if (!computeGradients) {*dmGrad = NULL; return(0);}
 78:   DMGetDMTS(dm, &dmts);
 79:   PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", &obj);
 80:   if (!obj) {
 81:     DM  dmGrad;
 82:     Vec faceGeometry, cellGeometry;

 84:     DMPlexTSGetGeometryFVM(dm, &faceGeometry, &cellGeometry, NULL);
 85:     DMPlexComputeGradientFVM(dm, fv, faceGeometry, cellGeometry, &dmGrad);
 86:     PetscObjectCompose((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject) dmGrad);
 87:     DMDestroy(&dmGrad);
 88:   }
 89:   PetscObjectQuery((PetscObject) dmts, "DMPlexTS_dmgrad_fvm", (PetscObject *) dmGrad);
 90:   return(0);
 91: }

 95: /*@
 96:   DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user

 98:   Input Parameters:
 99: + dm - The mesh
100: . t - The time
101: . locX  - Local solution
102: - user - The user context

104:   Output Parameter:
105: . F  - Global output vector

107:   Level: developer

109: .seealso: DMPlexComputeJacobianActionFEM()
110: @*/
111: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
112: {
113:   Vec            locF;
114:   PetscInt       cStart, cEnd, cEndInterior;

118:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
119:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
120:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
121:   DMGetLocalVector(dm, &locF);
122:   VecZeroEntries(locF);
123:   DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, NULL, locF, user);
124:   DMLocalToGlobalBegin(dm, locF, INSERT_VALUES, F);
125:   DMLocalToGlobalEnd(dm, locF, INSERT_VALUES, F);
126:   DMRestoreLocalVector(dm, &locF);
127:   return(0);
128: }

132: /*@
133:   DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user

135:   Input Parameters:
136: + dm - The mesh
137: . t - The time
138: . locX  - Local solution
139: . locX_t - Local solution time derivative, or NULL
140: - user - The user context

142:   Output Parameter:
143: . locF  - Local output vector

145:   Level: developer

147: .seealso: DMPlexComputeJacobianActionFEM()
148: @*/
149: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
150: {
151:   PetscInt       cStart, cEnd, cEndInterior;

155:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
156:   DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);
157:   cEnd = cEndInterior < 0 ? cEnd : cEndInterior;
158:   DMPlexComputeResidual_Internal(dm, cStart, cEnd, time, locX, locX_t, locF, user);
159:   return(0);
160: }

164: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u, PetscErrorCode (**exactFuncs)(PetscInt dim, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx), void **ctxs)
165: {
166:   DM             dm;
167:   SNES           snes;
168:   Vec            sol;
169:   PetscBool      check;

173:   PetscOptionsHasName(ts->hdr.prefix, "-dmts_check", &check);
174:   if (!check) return(0);
175:   VecDuplicate(u, &sol);
176:   TSSetSolution(ts, sol);
177:   TSGetDM(ts, &dm);
178:   TSSetUp(ts);
179:   TSGetSNES(ts, &snes);
180:   SNESSetSolution(snes, sol);
181:   DMSNESCheckFromOptions_Internal(snes, dm, u, sol, exactFuncs, ctxs);
182:   VecDestroy(&sol);
183:   return(0);
184: }