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;

 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:         DMCopyDMTS(dm, *plex);
 24:         DMCopyDMSNES(dm, *plex);
 25:         DMCopyAuxiliaryVec(dm, *plex);
 26:       }
 27:     } else {
 28:       PetscObjectReference((PetscObject) *plex);
 29:     }
 30:   }
 31:   return(0);
 32: }

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

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

 43:   Output Parameter:
 44: . F  - Global output vector

 46:   Level: developer

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

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

 77: /*@
 78:   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

 80:   Input Parameters:
 81: + dm - The mesh
 82: . t - The time
 83: . locX  - Local solution
 84: . locX_t - Local solution time derivative, or NULL
 85: - user - The user context

 87:   Level: developer

 89: .seealso: DMPlexComputeJacobianActionFEM()
 90: @*/
 91: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
 92: {
 93:   DM             plex;
 94:   Vec            faceGeometryFVM = NULL;
 95:   PetscInt       Nf, f;

 99:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
100:   DMGetNumFields(plex, &Nf);
101:   if (!locX_t) {
102:     /* This is the RHS part */
103:     for (f = 0; f < Nf; f++) {
104:       PetscObject  obj;
105:       PetscClassId id;

107:       DMGetField(plex, f, NULL, &obj);
108:       PetscObjectGetClassId(obj, &id);
109:       if (id == PETSCFV_CLASSID) {
110:         DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
111:         break;
112:       }
113:     }
114:   }
115:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
116:   DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
117:   DMDestroy(&plex);
118:   return(0);
119: }

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

124:   Input Parameters:
125: + dm - The mesh
126: . t - The time
127: . locX  - Local solution
128: . locX_t - Local solution time derivative, or NULL
129: - user - The user context

131:   Output Parameter:
132: . locF  - Local output vector

134:   Level: developer

136: .seealso: DMPlexComputeJacobianActionFEM()
137: @*/
138: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
139: {
140:   DM             plex;
141:   IS             allcellIS;
142:   PetscInt       Nds, s;

146:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
147:   DMPlexGetAllCells_Internal(plex, &allcellIS);
148:   DMGetNumDS(dm, &Nds);
149:   for (s = 0; s < Nds; ++s) {
150:     PetscDS          ds;
151:     IS               cellIS;
152:     PetscFormKey key;

154:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
155:     key.value = 0;
156:     key.field = 0;
157:     key.part  = 0;
158:     if (!key.label) {
159:       PetscObjectReference((PetscObject) allcellIS);
160:       cellIS = allcellIS;
161:     } else {
162:       IS pointIS;

164:       key.value = 1;
165:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
166:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
167:       ISDestroy(&pointIS);
168:     }
169:     DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user);
170:     ISDestroy(&cellIS);
171:   }
172:   ISDestroy(&allcellIS);
173:   DMDestroy(&plex);
174:   return(0);
175: }

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

180:   Input Parameters:
181: + dm - The mesh
182: . t - The time
183: . locX  - Local solution
184: . locX_t - Local solution time derivative, or NULL
185: . X_tshift - The multiplicative parameter for dF/du_t
186: - user - The user context

188:   Output Parameter:
189: . locF  - Local output vector

191:   Level: developer

193: .seealso: DMPlexComputeJacobianActionFEM()
194: @*/
195: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
196: {
197:   DM             plex;
198:   IS             allcellIS;
199:   PetscBool      hasJac, hasPrec;
200:   PetscInt       Nds, s;

204:   DMTSConvertPlex(dm, &plex, PETSC_TRUE);
205:   DMPlexGetAllCells_Internal(plex, &allcellIS);
206:   DMGetNumDS(dm, &Nds);
207:   for (s = 0; s < Nds; ++s) {
208:     PetscDS          ds;
209:     IS               cellIS;
210:     PetscFormKey key;

212:     DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
213:     key.value = 0;
214:     key.field = 0;
215:     key.part  = 0;
216:     if (!key.label) {
217:       PetscObjectReference((PetscObject) allcellIS);
218:       cellIS = allcellIS;
219:     } else {
220:       IS pointIS;

222:       key.value = 1;
223:       DMLabelGetStratumIS(key.label, key.value, &pointIS);
224:       ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
225:       ISDestroy(&pointIS);
226:     }
227:     if (!s) {
228:       PetscDSHasJacobian(ds, &hasJac);
229:       PetscDSHasJacobianPreconditioner(ds, &hasPrec);
230:       if (hasJac && hasPrec) {MatZeroEntries(Jac);}
231:       MatZeroEntries(JacP);
232:     }
233:     DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
234:     ISDestroy(&cellIS);
235:   }
236:   ISDestroy(&allcellIS);
237:   DMDestroy(&plex);
238:   return(0);
239: }

241: /*@C
242:   DMTSCheckResidual - Check the residual of the exact solution

244:   Input Parameters:
245: + ts  - the TS object
246: . dm  - the DM
247: . t   - the time
248: . u   - a DM vector
249: . u_t - a DM vector
250: - tol - A tolerance for the check, or -1 to print the results instead

252:   Output Parameters:
253: . residual - The residual norm of the exact solution, or NULL

255:   Level: developer

257: .seealso: DNTSCheckFromOptions(), DMTSCheckJacobian(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
258: @*/
259: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
260: {
261:   MPI_Comm       comm;
262:   Vec            r;
263:   PetscReal      res;

271:   PetscObjectGetComm((PetscObject) ts, &comm);
272:   DMComputeExactSolution(dm, t, u, u_t);
273:   VecDuplicate(u, &r);
274:   TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
275:   VecNorm(r, NORM_2, &res);
276:   if (tol >= 0.0) {
277:     if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
278:   } else if (residual) {
279:     *residual = res;
280:   } else {
281:     PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
282:     VecChop(r, 1.0e-10);
283:     PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) dm);
284:     PetscObjectSetName((PetscObject) r, "Initial Residual");
285:     PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
286:     VecViewFromOptions(r, NULL, "-vec_view");
287:     PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
288:   }
289:   VecDestroy(&r);
290:   return(0);
291: }

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

296:   Input Parameters:
297: + ts  - the TS object
298: . dm  - the DM
299: . t   - the time
300: . u   - a DM vector
301: . u_t - a DM vector
302: - tol - A tolerance for the check, or -1 to print the results instead

304:   Output Parameters:
305: + isLinear - Flag indicaing that the function looks linear, or NULL
306: - convRate - The rate of convergence of the linear model, or NULL

308:   Level: developer

310: .seealso: DNTSCheckFromOptions(), DMTSCheckResidual(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
311: @*/
312: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
313: {
314:   MPI_Comm       comm;
315:   PetscDS        ds;
316:   Mat            J, M;
317:   MatNullSpace   nullspace;
318:   PetscReal      dt, shift, slope, intercept;
319:   PetscBool      hasJac, hasPrec, isLin = PETSC_FALSE;

328:   PetscObjectGetComm((PetscObject) ts, &comm);
329:   DMComputeExactSolution(dm, t, u, u_t);
330:   /* Create and view matrices */
331:   TSGetTimeStep(ts, &dt);
332:   shift = 1.0/dt;
333:   DMCreateMatrix(dm, &J);
334:   DMGetDS(dm, &ds);
335:   PetscDSHasJacobian(ds, &hasJac);
336:   PetscDSHasJacobianPreconditioner(ds, &hasPrec);
337:   if (hasJac && hasPrec) {
338:     DMCreateMatrix(dm, &M);
339:     TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
340:     PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
341:     PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
342:     MatViewFromOptions(M, NULL, "-mat_view");
343:     MatDestroy(&M);
344:   } else {
345:     TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
346:   }
347:   PetscObjectSetName((PetscObject) J, "Jacobian");
348:   PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
349:   MatViewFromOptions(J, NULL, "-mat_view");
350:   /* Check nullspace */
351:   MatGetNullSpace(J, &nullspace);
352:   if (nullspace) {
353:     PetscBool isNull;
354:     MatNullSpaceTest(nullspace, J, &isNull);
355:     if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
356:   }
357:   /* Taylor test */
358:   {
359:     PetscRandom rand;
360:     Vec         du, uhat, uhat_t, r, rhat, df;
361:     PetscReal   h;
362:     PetscReal  *es, *hs, *errors;
363:     PetscReal   hMax = 1.0, hMin = 1e-6, hMult = 0.1;
364:     PetscInt    Nv, v;

366:     /* Choose a perturbation direction */
367:     PetscRandomCreate(comm, &rand);
368:     VecDuplicate(u, &du);
369:     VecSetRandom(du, rand);
370:     PetscRandomDestroy(&rand);
371:     VecDuplicate(u, &df);
372:     MatMult(J, du, df);
373:     /* Evaluate residual at u, F(u), save in vector r */
374:     VecDuplicate(u, &r);
375:     TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
376:     /* Look at the convergence of our Taylor approximation as we approach u */
377:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
378:     PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
379:     VecDuplicate(u, &uhat);
380:     VecDuplicate(u, &uhat_t);
381:     VecDuplicate(u, &rhat);
382:     for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
383:       VecWAXPY(uhat, h, du, u);
384:       VecWAXPY(uhat_t, h*shift, du, u_t);
385:       /* 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 */
386:       TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
387:       VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
388:       VecNorm(rhat, NORM_2, &errors[Nv]);

390:       es[Nv] = PetscLog10Real(errors[Nv]);
391:       hs[Nv] = PetscLog10Real(h);
392:     }
393:     VecDestroy(&uhat);
394:     VecDestroy(&uhat_t);
395:     VecDestroy(&rhat);
396:     VecDestroy(&df);
397:     VecDestroy(&r);
398:     VecDestroy(&du);
399:     for (v = 0; v < Nv; ++v) {
400:       if ((tol >= 0) && (errors[v] > tol)) break;
401:       else if (errors[v] > PETSC_SMALL)    break;
402:     }
403:     if (v == Nv) isLin = PETSC_TRUE;
404:     PetscLinearRegression(Nv, hs, es, &slope, &intercept);
405:     PetscFree3(es, hs, errors);
406:     /* Slope should be about 2 */
407:     if (tol >= 0) {
408:       if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
409:     } else if (isLinear || convRate) {
410:       if (isLinear) *isLinear = isLin;
411:       if (convRate) *convRate = slope;
412:     } else {
413:       if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
414:       else        {PetscPrintf(comm, "Function appears to be linear\n");}
415:     }
416:   }
417:   MatDestroy(&J);
418:   return(0);
419: }

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

424:   Input Parameters:
425: + ts - the TS object
426: - u  - representative TS vector

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

430:   Level: developer
431: @*/
432: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
433: {
434:   DM             dm;
435:   SNES           snes;
436:   Vec            sol, u_t;
437:   PetscReal      t;
438:   PetscBool      check;

442:   PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
443:   if (!check) return(0);
444:   VecDuplicate(u, &sol);
445:   VecCopy(u, sol);
446:   TSSetSolution(ts, u);
447:   TSGetDM(ts, &dm);
448:   TSSetUp(ts);
449:   TSGetSNES(ts, &snes);
450:   SNESSetSolution(snes, u);

452:   TSGetTime(ts, &t);
453:   DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
454:   DMGetGlobalVector(dm, &u_t);
455:   DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
456:   DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
457:   DMRestoreGlobalVector(dm, &u_t);

459:   VecDestroy(&sol);
460:   return(0);
461: }