Actual source code: dmplexts.c

petsc-master 2019-06-22
Report Typos and Errors
  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:   DMSNESCheck_Internal(snes, dm, sol, exactFuncs, ctxs);
285:   VecDestroy(&sol);
286:   return(0);
287: }