Actual source code: dmplexts.c

petsc-master 2020-08-06
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: }

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

 44:   Input Parameters:
 45: + dm - The mesh
 46: . t - The time
 47: . locX  - Local solution
 48: - user - The user context

 50:   Output Parameter:
 51: . F  - Global output vector

 53:   Level: developer

 55: .seealso: DMPlexComputeJacobianActionFEM()
 56: @*/
 57: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
 58: {
 59:   Vec            locF;
 60:   IS             cellIS;
 61:   DM             plex;
 62:   PetscInt       depth;

 66:   DMTSConvertPlex(dm,&plex,PETSC_TRUE);
 67:   DMPlexGetDepth(plex, &depth);
 68:   DMGetStratumIS(plex, "dim", depth, &cellIS);
 69:   if (!cellIS) {
 70:     DMGetStratumIS(plex, "depth", depth, &cellIS);
 71:   }
 72:   DMGetLocalVector(plex, &locF);
 73:   VecZeroEntries(locF);
 74:   DMPlexComputeResidual_Internal(plex, cellIS, time, locX, NULL, time, locF, user);
 75:   DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
 76:   DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
 77:   DMRestoreLocalVector(plex, &locF);
 78:   ISDestroy(&cellIS);
 79:   DMDestroy(&plex);
 80:   return(0);
 81: }

 83: /*@
 84:   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

 86:   Input Parameters:
 87: + dm - The mesh
 88: . t - The time
 89: . locX  - Local solution
 90: . locX_t - Local solution time derivative, or NULL
 91: - user - The user context

 93:   Level: developer

 95: .seealso: DMPlexComputeJacobianActionFEM()
 96: @*/
 97: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
 98: {
 99:   DM             plex;
100:   Vec            faceGeometryFVM = NULL;
101:   PetscInt       Nf, f;

105:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
106:   DMGetNumFields(plex, &Nf);
107:   if (!locX_t) {
108:     /* This is the RHS part */
109:     for (f = 0; f < Nf; f++) {
110:       PetscObject  obj;
111:       PetscClassId id;

113:       DMGetField(plex, f, NULL, &obj);
114:       PetscObjectGetClassId(obj, &id);
115:       if (id == PETSCFV_CLASSID) {
116:         DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
117:         break;
118:       }
119:     }
120:   }
121:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
122:   DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
123:   DMDestroy(&plex);
124:   return(0);
125: }

127: /*@
128:   DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X 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:   Output Parameter:
138: . locF  - Local output vector

140:   Level: developer

142: .seealso: DMPlexComputeJacobianActionFEM()
143: @*/
144: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
145: {
146:   DM             plex;
147:   IS             cellIS;
148:   PetscInt       depth;

152:   DMTSConvertPlex(dm,&plex,PETSC_TRUE);
153:   DMPlexGetDepth(plex, &depth);
154:   DMGetStratumIS(plex, "dim", depth, &cellIS);
155:   if (!cellIS) {
156:     DMGetStratumIS(plex, "depth", depth, &cellIS);
157:   }
158:   DMPlexComputeResidual_Internal(plex, cellIS, time, locX, locX_t, time, locF, user);
159:   ISDestroy(&cellIS);
160:   DMDestroy(&plex);
161:   return(0);
162: }

164: /*@
165:   DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user

167:   Input Parameters:
168: + dm - The mesh
169: . t - The time
170: . locX  - Local solution
171: . locX_t - Local solution time derivative, or NULL
172: . X_tshift - The multiplicative parameter for dF/du_t
173: - user - The user context

175:   Output Parameter:
176: . locF  - Local output vector

178:   Level: developer

180: .seealso: DMPlexComputeJacobianActionFEM()
181: @*/
182: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
183: {
184:   DM             plex;
185:   PetscDS        prob;
186:   PetscBool      hasJac, hasPrec;
187:   IS             cellIS;
188:   PetscInt       depth;

192:   DMTSConvertPlex(dm,&plex,PETSC_TRUE);
193:   DMGetDS(dm, &prob);
194:   PetscDSHasJacobian(prob, &hasJac);
195:   PetscDSHasJacobianPreconditioner(prob, &hasPrec);
196:   if (hasJac && hasPrec) {MatZeroEntries(Jac);}
197:   MatZeroEntries(JacP);
198:   DMPlexGetDepth(plex,&depth);
199:   DMGetStratumIS(plex, "dim", depth, &cellIS);
200:   if (!cellIS) {DMGetStratumIS(plex, "depth", depth, &cellIS);}
201:   DMPlexComputeJacobian_Internal(plex, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
202:   ISDestroy(&cellIS);
203:   DMDestroy(&plex);
204:   return(0);
205: }

207: /*@C
208:   DMTSCheckResidual - Check the residual of the exact solution

210:   Input Parameters:
211: + ts  - the TS object
212: . dm  - the DM
213: . t   - the time
214: . u   - a DM vector
215: . u_t - a DM vector
216: - tol - A tolerance for the check, or -1 to print the results instead

218:   Output Parameters:
219: . residual - The residual norm of the exact solution, or NULL

221:   Level: developer

223: .seealso: DNTSCheckFromOptions(), DMTSCheckJacobian(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
224: @*/
225: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
226: {
227:   MPI_Comm       comm;
228:   Vec            r;
229:   PetscReal      res;

237:   PetscObjectGetComm((PetscObject) ts, &comm);
238:   DMComputeExactSolution(dm, t, u, u_t);
239:   VecDuplicate(u, &r);
240:   TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
241:   VecNorm(r, NORM_2, &res);
242:   if (tol >= 0.0) {
243:     if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
244:   } else if (residual) {
245:     *residual = res;
246:   } else {
247:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
248:     VecChop(r, 1.0e-10);
249:     PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) dm);
250:     PetscObjectSetName((PetscObject) r, "Initial Residual");
251:     PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
252:     VecViewFromOptions(r, NULL, "-vec_view");
253:     PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
254:   }
255:   VecDestroy(&r);
256:   return(0);
257: }

259: /*@C
260:   DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test

262:   Input Parameters:
263: + ts  - the TS object
264: . dm  - the DM
265: . t   - the time
266: . u   - a DM vector
267: . u_t - a DM vector
268: - tol - A tolerance for the check, or -1 to print the results instead

270:   Output Parameters:
271: + isLinear - Flag indicaing that the function looks linear, or NULL
272: - convRate - The rate of convergence of the linear model, or NULL

274:   Level: developer

276: .seealso: DNTSCheckFromOptions(), DMTSCheckResidual(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
277: @*/
278: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
279: {
280:   MPI_Comm       comm;
281:   PetscDS        ds;
282:   Mat            J, M;
283:   MatNullSpace   nullspace;
284:   PetscReal      dt, shift, slope, intercept;
285:   PetscBool      hasJac, hasPrec, isLin = PETSC_FALSE;

294:   PetscObjectGetComm((PetscObject) ts, &comm);
295:   DMComputeExactSolution(dm, t, u, u_t);
296:   /* Create and view matrices */
297:   TSGetTimeStep(ts, &dt);
298:   shift = 1.0/dt;
299:   DMCreateMatrix(dm, &J);
300:   DMGetDS(dm, &ds);
301:   PetscDSHasJacobian(ds, &hasJac);
302:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
303:   if (hasJac && hasPrec) {
304:     DMCreateMatrix(dm, &M);
305:     TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
306:     PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
307:     PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
308:     MatViewFromOptions(M, NULL, "-mat_view");
309:     MatDestroy(&M);
310:   } else {
311:     TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
312:   }
313:   PetscObjectSetName((PetscObject) J, "Jacobian");
314:   PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
315:   MatViewFromOptions(J, NULL, "-mat_view");
316:   /* Check nullspace */
317:   MatGetNullSpace(J, &nullspace);
318:   if (nullspace) {
319:     PetscBool isNull;
320:     MatNullSpaceTest(nullspace, J, &isNull);
321:     if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
322:   }
323:   MatNullSpaceDestroy(&nullspace);
324:   /* Taylor test */
325:   {
326:     PetscRandom rand;
327:     Vec         du, uhat, uhat_t, r, rhat, df;
328:     PetscReal   h;
329:     PetscReal  *es, *hs, *errors;
330:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
331:     PetscInt    Nv, v;

333:     /* Choose a perturbation direction */
334:     PetscRandomCreate(comm, &rand);
335:     VecDuplicate(u, &du);
336:     VecSetRandom(du, rand); 
337:     PetscRandomDestroy(&rand);
338:     VecDuplicate(u, &df);
339:     MatMult(J, du, df);
340:     /* Evaluate residual at u, F(u), save in vector r */
341:     VecDuplicate(u, &r);
342:     TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
343:     /* Look at the convergence of our Taylor approximation as we approach u */
344:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
345:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
346:     VecDuplicate(u, &uhat);
347:     VecDuplicate(u, &uhat_t);
348:     VecDuplicate(u, &rhat);
349:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
350:       VecWAXPY(uhat, h, du, u);
351:       VecWAXPY(uhat_t, h*shift, du, u_t);
352:       /* F(\hat u, \hat u_t) \approx F(u, u_t) + J(u, u_t) (uhat - u) + J_t(u, u_t) (uhat_t - u_t) = F(u) + h * J(u) du + h * shift * J_t(u) du = F(u) + h F' du */
353:       TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
354:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
355:       VecNorm(rhat, NORM_2, &errors[Nv]);

357:       es[Nv] = PetscLog10Real(errors[Nv]);
358:       hs[Nv] = PetscLog10Real(h);
359:     }
360:     VecDestroy(&uhat);
361:     VecDestroy(&uhat_t);
362:     VecDestroy(&rhat);
363:     VecDestroy(&df);
364:     VecDestroy(&r);
365:     VecDestroy(&du);
366:     for (v = 0; v < Nv; ++v) {
367:       if ((tol >= 0) && (errors[v] > tol)) break;
368:       else if (errors[v] > PETSC_SMALL)    break;
369:     }
370:     if (v == Nv) isLin = PETSC_TRUE;
371:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
372:     PetscFree3(es, hs, errors);
373:     /* Slope should be about 2 */
374:     if (tol >= 0) {
375:       if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
376:     } else if (isLinear || convRate) {
377:       if (isLinear) *isLinear = isLin;
378:       if (convRate) *convRate = slope;
379:     } else {
380:       if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
381:       else        {PetscPrintf(comm, "Function appears to be linear\n");}
382:     }
383:   }
384:   MatDestroy(&J);
385:   return(0);
386: }

388: /*@C
389:   DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information

391:   Input Parameters:
392: + ts - the TS object
393: - u  - representative TS vector

395:   Note: The user must call PetscDSSetExactSolution() beforehand

397:   Level: developer
398: @*/
399: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
400: {
401:   DM             dm;
402:   SNES           snes;
403:   Vec            sol, u_t;
404:   PetscReal      t;
405:   PetscBool      check;

409:   PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
410:   if (!check) return(0);
411:   VecDuplicate(u, &sol);
412:   VecCopy(u, sol);
413:   TSSetSolution(ts, u);
414:   TSGetDM(ts, &dm);
415:   TSSetUp(ts);
416:   TSGetSNES(ts, &snes);
417:   SNESSetSolution(snes, u);

419:   TSGetTime(ts, &t);
420:   DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
421:   DMGetGlobalVector(dm, &u_t);
422:   DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
423:   DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
424:   DMRestoreGlobalVector(dm, &u_t);

426:   VecDestroy(&sol);
427:   return(0);
428: }