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: 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: 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: 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: 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: 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: DNTSCheckFromOptions(), 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: 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:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
416:     VecDuplicate(u, &uhat);
417:     VecDuplicate(u, &uhat_t);
418:     VecDuplicate(u, &rhat);
419:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
420:       VecWAXPY(uhat, h, du, u);
421:       VecWAXPY(uhat_t, h*shift, du, u_t);
422:       /* 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 */
423:       TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
424:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
425:       VecNorm(rhat, NORM_2, &errors[Nv]);

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

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

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

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

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

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

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

494:   VecDestroy(&sol);
495:   return 0;
496: }