Actual source code: dmplexts.c

  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;

 11:   PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex);
 12:   if (isPlex) {
 13:     *plex = dm;
 14:     PetscObjectReference((PetscObject)dm);
 15:   } else {
 16:     PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex);
 17:     if (!*plex) {
 18:       DMConvert(dm, DMPLEX, plex);
 19:       PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex);
 20:       if (copy) {
 21:         DMCopyDMTS(dm, *plex);
 22:         DMCopyDMSNES(dm, *plex);
 23:         DMCopyAuxiliaryVec(dm, *plex);
 24:       }
 25:     } else {
 26:       PetscObjectReference((PetscObject)*plex);
 27:     }
 28:   }
 29:   return 0;
 30: }

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

 35:   Input Parameters:
 36: + dm - The mesh
 37: . t - The time
 38: . locX  - Local solution
 39: - user - The user context

 41:   Output Parameter:
 42: . F  - Global output vector

 44:   Level: developer

 46: .seealso: [](chapter_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
 47: @*/
 48: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
 49: {
 50:   Vec          locF;
 51:   IS           cellIS;
 52:   DM           plex;
 53:   PetscInt     depth;
 54:   PetscFormKey key = {NULL, 0, 0, 0};

 56:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
 57:   DMPlexGetDepth(plex, &depth);
 58:   DMGetStratumIS(plex, "dim", depth, &cellIS);
 59:   if (!cellIS) DMGetStratumIS(plex, "depth", depth, &cellIS);
 60:   DMGetLocalVector(plex, &locF);
 61:   VecZeroEntries(locF);
 62:   DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locF, user);
 63:   DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
 64:   DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
 65:   DMRestoreLocalVector(plex, &locF);
 66:   ISDestroy(&cellIS);
 67:   DMDestroy(&plex);
 68:   return 0;
 69: }

 71: /*@
 72:   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

 74:   Input Parameters:
 75: + dm - The mesh
 76: . t - The time
 77: . locX  - Local solution
 78: . locX_t - Local solution time derivative, or NULL
 79: - user - The user context

 81:   Level: developer

 83: .seealso: [](chapter_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
 84: @*/
 85: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
 86: {
 87:   DM       plex;
 88:   Vec      faceGeometryFVM = NULL;
 89:   PetscInt Nf, f;

 91:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
 92:   DMGetNumFields(plex, &Nf);
 93:   if (!locX_t) {
 94:     /* This is the RHS part */
 95:     for (f = 0; f < Nf; f++) {
 96:       PetscObject  obj;
 97:       PetscClassId id;

 99:       DMGetField(plex, f, NULL, &obj);
100:       PetscObjectGetClassId(obj, &id);
101:       if (id == PETSCFV_CLASSID) {
102:         DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
103:         break;
104:       }
105:     }
106:   }
107:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
108:   DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
109:   DMDestroy(&plex);
110:   return 0;
111: }

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

116:   Input Parameters:
117: + dm - The mesh
118: . t - The time
119: . locX  - Local solution
120: . locX_t - Local solution time derivative, or NULL
121: - user - The user context

123:   Output Parameter:
124: . locF  - Local output vector

126:   Level: developer

128: .seealso: [](chapter_ts), `DMPLEX`, `TS`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
129: @*/
130: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
131: {
132:   DM       plex;
133:   IS       allcellIS;
134:   PetscInt Nds, s;

136:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
137:   DMPlexGetAllCells_Internal(plex, &allcellIS);
138:   DMGetNumDS(dm, &Nds);
139:   for (s = 0; s < Nds; ++s) {
140:     PetscDS      ds;
141:     IS           cellIS;
142:     PetscFormKey key;

144:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
145:     key.value = 0;
146:     key.field = 0;
147:     key.part  = 0;
148:     if (!key.label) {
149:       PetscObjectReference((PetscObject)allcellIS);
150:       cellIS = allcellIS;
151:     } else {
152:       IS pointIS;

154:       key.value = 1;
155:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
156:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
157:       ISDestroy(&pointIS);
158:     }
159:     DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user);
160:     ISDestroy(&cellIS);
161:   }
162:   ISDestroy(&allcellIS);
163:   DMDestroy(&plex);
164:   return 0;
165: }

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

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

178:   Output Parameter:
179: . locF  - Local output vector

181:   Level: developer

183: .seealso: [](chapter_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
184: @*/
185: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
186: {
187:   DM        plex;
188:   IS        allcellIS;
189:   PetscBool hasJac, hasPrec;
190:   PetscInt  Nds, s;

192:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
193:   DMPlexGetAllCells_Internal(plex, &allcellIS);
194:   DMGetNumDS(dm, &Nds);
195:   for (s = 0; s < Nds; ++s) {
196:     PetscDS      ds;
197:     IS           cellIS;
198:     PetscFormKey key;

200:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
201:     key.value = 0;
202:     key.field = 0;
203:     key.part  = 0;
204:     if (!key.label) {
205:       PetscObjectReference((PetscObject)allcellIS);
206:       cellIS = allcellIS;
207:     } else {
208:       IS pointIS;

210:       key.value = 1;
211:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
212:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
213:       ISDestroy(&pointIS);
214:     }
215:     if (!s) {
216:       PetscDSHasJacobian(ds, &hasJac);
217:       PetscDSHasJacobianPreconditioner(ds, &hasPrec);
218:       if (hasJac && hasPrec) MatZeroEntries(Jac);
219:       MatZeroEntries(JacP);
220:     }
221:     DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
222:     ISDestroy(&cellIS);
223:   }
224:   ISDestroy(&allcellIS);
225:   DMDestroy(&plex);
226:   return 0;
227: }

229: /*@
230:   DMPlexTSComputeRHSFunctionFEM - Form the local residual G from the local input X using pointwise functions specified by the user

232:   Input Parameters:
233: + dm - The mesh
234: . t - The time
235: . locX  - Local solution
236: - user - The user context

238:   Output Parameter:
239: . locG  - Local output vector

241:   Level: developer

243: .seealso: [](chapter_ts),  `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeIJacobianFEM()`
244: @*/
245: PetscErrorCode DMPlexTSComputeRHSFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locG, void *user)
246: {
247:   DM       plex;
248:   IS       allcellIS;
249:   PetscInt Nds, s;

251:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
252:   DMPlexGetAllCells_Internal(plex, &allcellIS);
253:   DMGetNumDS(dm, &Nds);
254:   for (s = 0; s < Nds; ++s) {
255:     PetscDS      ds;
256:     IS           cellIS;
257:     PetscFormKey key;

259:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
260:     key.value = 0;
261:     key.field = 0;
262:     key.part  = 100;
263:     if (!key.label) {
264:       PetscObjectReference((PetscObject)allcellIS);
265:       cellIS = allcellIS;
266:     } else {
267:       IS pointIS;

269:       key.value = 1;
270:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
271:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
272:       ISDestroy(&pointIS);
273:     }
274:     DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locG, user);
275:     ISDestroy(&cellIS);
276:   }
277:   ISDestroy(&allcellIS);
278:   DMDestroy(&plex);
279:   return 0;
280: }

282: /*@C
283:   DMTSCheckResidual - Check the residual of the exact solution

285:   Input Parameters:
286: + ts  - the `TS` object
287: . dm  - the `DM`
288: . t   - the time
289: . u   - a `DM` vector
290: . u_t - a `DM` vector
291: - tol - A tolerance for the check, or -1 to print the results instead

293:   Output Parameters:
294: . residual - The residual norm of the exact solution, or NULL

296:   Level: developer

298: .seealso: [](chapter_ts), `DM`, `DMTSCheckFromOptions()`, `DMTSCheckJacobian()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
299: @*/
300: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
301: {
302:   MPI_Comm  comm;
303:   Vec       r;
304:   PetscReal res;

310:   PetscObjectGetComm((PetscObject)ts, &comm);
311:   DMComputeExactSolution(dm, t, u, u_t);
312:   VecDuplicate(u, &r);
313:   TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
314:   VecNorm(r, NORM_2, &res);
315:   if (tol >= 0.0) {
317:   } else if (residual) {
318:     *residual = res;
319:   } else {
320:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
321:     VecChop(r, 1.0e-10);
322:     PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)dm);
323:     PetscObjectSetName((PetscObject)r, "Initial Residual");
324:     PetscObjectSetOptionsPrefix((PetscObject)r, "res_");
325:     VecViewFromOptions(r, NULL, "-vec_view");
326:     PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL);
327:   }
328:   VecDestroy(&r);
329:   return 0;
330: }

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

335:   Input Parameters:
336: + ts  - the TS object
337: . dm  - the DM
338: . t   - the time
339: . u   - a DM vector
340: . u_t - a DM vector
341: - tol - A tolerance for the check, or -1 to print the results instead

343:   Output Parameters:
344: + isLinear - Flag indicaing that the function looks linear, or NULL
345: - convRate - The rate of convergence of the linear model, or NULL

347:   Level: developer

349: .seealso: [](chapter_ts), `DNTSCheckFromOptions()`, `DMTSCheckResidual()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
350: @*/
351: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
352: {
353:   MPI_Comm     comm;
354:   PetscDS      ds;
355:   Mat          J, M;
356:   MatNullSpace nullspace;
357:   PetscReal    dt, shift, slope, intercept;
358:   PetscBool    hasJac, hasPrec, isLin = PETSC_FALSE;

365:   PetscObjectGetComm((PetscObject)ts, &comm);
366:   DMComputeExactSolution(dm, t, u, u_t);
367:   /* Create and view matrices */
368:   TSGetTimeStep(ts, &dt);
369:   shift = 1.0 / dt;
370:   DMCreateMatrix(dm, &J);
371:   DMGetDS(dm, &ds);
372:   PetscDSHasJacobian(ds, &hasJac);
373:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
374:   if (hasJac && hasPrec) {
375:     DMCreateMatrix(dm, &M);
376:     TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
377:     PetscObjectSetName((PetscObject)M, "Preconditioning Matrix");
378:     PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_");
379:     MatViewFromOptions(M, NULL, "-mat_view");
380:     MatDestroy(&M);
381:   } else {
382:     TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
383:   }
384:   PetscObjectSetName((PetscObject)J, "Jacobian");
385:   PetscObjectSetOptionsPrefix((PetscObject)J, "jac_");
386:   MatViewFromOptions(J, NULL, "-mat_view");
387:   /* Check nullspace */
388:   MatGetNullSpace(J, &nullspace);
389:   if (nullspace) {
390:     PetscBool isNull;
391:     MatNullSpaceTest(nullspace, J, &isNull);
393:   }
394:   /* Taylor test */
395:   {
396:     PetscRandom rand;
397:     Vec         du, uhat, uhat_t, r, rhat, df;
398:     PetscReal   h;
399:     PetscReal  *es, *hs, *errors;
400:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
401:     PetscInt    Nv, v;

403:     /* Choose a perturbation direction */
404:     PetscRandomCreate(comm, &rand);
405:     VecDuplicate(u, &du);
406:     VecSetRandom(du, rand);
407:     PetscRandomDestroy(&rand);
408:     VecDuplicate(u, &df);
409:     MatMult(J, du, df);
410:     /* Evaluate residual at u, F(u), save in vector r */
411:     VecDuplicate(u, &r);
412:     TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
413:     /* Look at the convergence of our Taylor approximation as we approach u */
414:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
415:       ;
416:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
417:     VecDuplicate(u, &uhat);
418:     VecDuplicate(u, &uhat_t);
419:     VecDuplicate(u, &rhat);
420:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
421:       VecWAXPY(uhat, h, du, u);
422:       VecWAXPY(uhat_t, h * shift, du, u_t);
423:       /* 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 */
424:       TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
425:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
426:       VecNorm(rhat, NORM_2, &errors[Nv]);

428:       es[Nv] = PetscLog10Real(errors[Nv]);
429:       hs[Nv] = PetscLog10Real(h);
430:     }
431:     VecDestroy(&uhat);
432:     VecDestroy(&uhat_t);
433:     VecDestroy(&rhat);
434:     VecDestroy(&df);
435:     VecDestroy(&r);
436:     VecDestroy(&du);
437:     for (v = 0; v < Nv; ++v) {
438:       if ((tol >= 0) && (errors[v] > tol)) break;
439:       else if (errors[v] > PETSC_SMALL) break;
440:     }
441:     if (v == Nv) isLin = PETSC_TRUE;
442:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
443:     PetscFree3(es, hs, errors);
444:     /* Slope should be about 2 */
445:     if (tol >= 0) {
447:     } else if (isLinear || convRate) {
448:       if (isLinear) *isLinear = isLin;
449:       if (convRate) *convRate = slope;
450:     } else {
451:       if (!isLin) PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope);
452:       else PetscPrintf(comm, "Function appears to be linear\n");
453:     }
454:   }
455:   MatDestroy(&J);
456:   return 0;
457: }

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

462:   Input Parameters:
463: + ts - the `TS` object
464: - u  - representative `TS` vector

466:   Note:
467:   The user must call `PetscDSSetExactSolution()` beforehand

469:   Level: developer
470: @*/
471: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
472: {
473:   DM        dm;
474:   SNES      snes;
475:   Vec       sol, u_t;
476:   PetscReal t;
477:   PetscBool check;

479:   PetscOptionsHasName(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-dmts_check", &check);
480:   if (!check) return 0;
481:   VecDuplicate(u, &sol);
482:   VecCopy(u, sol);
483:   TSSetSolution(ts, u);
484:   TSGetDM(ts, &dm);
485:   TSSetUp(ts);
486:   TSGetSNES(ts, &snes);
487:   SNESSetSolution(snes, u);

489:   TSGetTime(ts, &t);
490:   DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
491:   DMGetGlobalVector(dm, &u_t);
492:   DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
493:   DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
494:   DMRestoreGlobalVector(dm, &u_t);

496:   VecDestroy(&sol);
497:   return 0;
498: }