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: PetscFunctionBegin;
12: PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex));
13: if (isPlex) {
14: *plex = dm;
15: PetscCall(PetscObjectReference((PetscObject)dm));
16: } else {
17: PetscCall(PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex));
18: if (!*plex) {
19: PetscCall(DMConvert(dm, DMPLEX, plex));
20: PetscCall(PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex));
21: } else {
22: PetscCall(PetscObjectReference((PetscObject)*plex));
23: }
24: if (copy) {
25: PetscCall(DMCopyDMTS(dm, *plex));
26: PetscCall(DMCopyDMSNES(dm, *plex));
27: PetscCall(DMCopyAuxiliaryVec(dm, *plex));
28: }
29: }
30: PetscFunctionReturn(PETSC_SUCCESS);
31: }
33: PetscErrorCode DMPlexTSComputeRHSFunctionFVMCEED(DM dm, PetscReal time, Vec locX, Vec F, void *user)
34: {
35: #ifdef PETSC_HAVE_LIBCEED
36: PetscFV fv;
37: Vec locF;
38: Ceed ceed;
39: DMCeed sd = dm->dmceed;
40: CeedVector clocX, clocF;
41: #endif
43: #ifdef PETSC_HAVE_LIBCEED
44: PetscFunctionBegin;
45: PetscCall(DMGetCeed(dm, &ceed));
46: PetscCheck(sd, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "This DM has no CEED data. Call DMCeedCreate() before computing the residual.");
47: if (time == 0.) PetscCall(DMCeedComputeGeometry(dm, sd));
48: PetscCall(DMGetField(dm, 0, NULL, (PetscObject *)&fv));
49: PetscCall(DMPlexInsertBoundaryValuesFVM(dm, fv, locX, time, NULL));
50: PetscCall(DMGetLocalVector(dm, &locF));
51: PetscCall(VecZeroEntries(locF));
52: PetscCall(VecGetCeedVectorRead(locX, ceed, &clocX));
53: PetscCall(VecGetCeedVector(locF, ceed, &clocF));
54: PetscCallCEED(CeedOperatorApplyAdd(sd->op, clocX, clocF, CEED_REQUEST_IMMEDIATE));
55: PetscCall(VecRestoreCeedVectorRead(locX, &clocX));
56: PetscCall(VecRestoreCeedVector(locF, &clocF));
57: PetscCall(DMLocalToGlobalBegin(dm, locF, ADD_VALUES, F));
58: PetscCall(DMLocalToGlobalEnd(dm, locF, ADD_VALUES, F));
59: PetscCall(DMRestoreLocalVector(dm, &locF));
60: PetscCall(VecViewFromOptions(F, NULL, "-fv_rhs_view"));
61: PetscFunctionReturn(PETSC_SUCCESS);
62: #else
63: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This requires libCEED. Reconfigure using --download-libceed");
64: #endif
65: }
67: /*@
68: DMPlexTSComputeRHSFunctionFVM - Form the forcing `F` from the local input `locX` using pointwise functions specified by the user
70: Input Parameters:
71: + dm - The mesh
72: . time - The time
73: . locX - Local solution
74: - user - The user context
76: Output Parameter:
77: . F - Global output vector
79: Level: developer
81: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
82: @*/
83: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
84: {
85: Vec locF;
86: IS cellIS;
87: DM plex;
88: PetscInt depth;
89: PetscFormKey key = {NULL, 0, 0, 0};
91: PetscFunctionBegin;
92: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
93: PetscCall(DMPlexGetDepth(plex, &depth));
94: PetscCall(DMGetStratumIS(plex, "dim", depth, &cellIS));
95: if (!cellIS) PetscCall(DMGetStratumIS(plex, "depth", depth, &cellIS));
96: PetscCall(DMGetLocalVector(plex, &locF));
97: PetscCall(VecZeroEntries(locF));
98: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locF, user));
99: PetscCall(DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F));
100: PetscCall(DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F));
101: PetscCall(DMRestoreLocalVector(plex, &locF));
102: PetscCall(ISDestroy(&cellIS));
103: PetscCall(DMDestroy(&plex));
104: PetscCall(VecViewFromOptions(F, NULL, "-fv_rhs_view"));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: /*@
109: DMPlexTSComputeBoundary - Insert the essential boundary values into the local input `locX` and/or its time derivative `locX_t` using pointwise functions specified by the user
111: Input Parameters:
112: + dm - The mesh
113: . time - The time
114: . locX - Local solution
115: . locX_t - Local solution time derivative, or `NULL`
116: - user - The user context
118: Level: developer
120: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
121: @*/
122: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
123: {
124: DM plex;
125: Vec faceGeometryFVM = NULL;
126: PetscInt Nf, f;
128: PetscFunctionBegin;
129: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
130: PetscCall(DMGetNumFields(plex, &Nf));
131: if (!locX_t) {
132: /* This is the RHS part */
133: for (f = 0; f < Nf; f++) {
134: PetscObject obj;
135: PetscClassId id;
137: PetscCall(DMGetField(plex, f, NULL, &obj));
138: PetscCall(PetscObjectGetClassId(obj, &id));
139: if (id == PETSCFV_CLASSID) {
140: PetscCall(DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL));
141: break;
142: }
143: }
144: }
145: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL));
146: PetscCall(DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL));
147: PetscCall(DMDestroy(&plex));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: /*@
152: DMPlexTSComputeIFunctionFEM - Form the local residual `locF` from the local input `locX` using pointwise functions specified by the user
154: Input Parameters:
155: + dm - The mesh
156: . time - The time
157: . locX - Local solution
158: . locX_t - Local solution time derivative, or `NULL`
159: - user - The user context
161: Output Parameter:
162: . locF - Local output vector
164: Level: developer
166: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexTSComputeRHSFunctionFEM()`
167: @*/
168: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
169: {
170: DM plex;
171: IS allcellIS;
172: PetscInt Nds, s;
174: PetscFunctionBegin;
175: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
176: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
177: PetscCall(DMGetNumDS(dm, &Nds));
178: for (s = 0; s < Nds; ++s) {
179: PetscDS ds;
180: IS cellIS;
181: PetscFormKey key;
183: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
184: key.value = 0;
185: key.field = 0;
186: key.part = 0;
187: if (!key.label) {
188: PetscCall(PetscObjectReference((PetscObject)allcellIS));
189: cellIS = allcellIS;
190: } else {
191: IS pointIS;
193: key.value = 1;
194: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
195: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
196: PetscCall(ISDestroy(&pointIS));
197: }
198: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user));
199: PetscCall(ISDestroy(&cellIS));
200: }
201: PetscCall(ISDestroy(&allcellIS));
202: PetscCall(DMDestroy(&plex));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*@
207: DMPlexTSComputeIJacobianFEM - Form the Jacobian `Jac` from the local input `locX` using pointwise functions specified by the user
209: Input Parameters:
210: + dm - The mesh
211: . time - The time
212: . locX - Local solution
213: . locX_t - Local solution time derivative, or `NULL`
214: . X_tShift - The multiplicative parameter for dF/du_t
215: - user - The user context
217: Output Parameters:
218: + Jac - the Jacobian
219: - JacP - an additional approximation for the Jacobian to be used to compute the preconditioner (often is `Jac`)
221: Level: developer
223: .seealso: [](ch_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
224: @*/
225: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
226: {
227: DM plex;
228: IS allcellIS;
229: PetscBool hasJac, hasPrec;
230: PetscInt Nds, s;
232: PetscFunctionBegin;
233: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
234: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
235: PetscCall(DMGetNumDS(dm, &Nds));
236: for (s = 0; s < Nds; ++s) {
237: PetscDS ds;
238: IS cellIS;
239: PetscFormKey key;
241: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
242: key.value = 0;
243: key.field = 0;
244: key.part = 0;
245: if (!key.label) {
246: PetscCall(PetscObjectReference((PetscObject)allcellIS));
247: cellIS = allcellIS;
248: } else {
249: IS pointIS;
251: key.value = 1;
252: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
253: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
254: PetscCall(ISDestroy(&pointIS));
255: }
256: if (!s) {
257: PetscCall(PetscDSHasJacobian(ds, &hasJac));
258: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
259: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
260: PetscCall(MatZeroEntries(JacP));
261: }
262: PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user));
263: PetscCall(ISDestroy(&cellIS));
264: }
265: PetscCall(ISDestroy(&allcellIS));
266: PetscCall(DMDestroy(&plex));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: /*@
271: DMPlexTSComputeRHSFunctionFEM - Form the local residual `locG` from the local input `locX` using pointwise functions specified by the user
273: Input Parameters:
274: + dm - The mesh
275: . time - The time
276: . locX - Local solution
277: - user - The user context
279: Output Parameter:
280: . locG - Local output vector
282: Level: developer
284: .seealso: [](ch_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeIJacobianFEM()`
285: @*/
286: PetscErrorCode DMPlexTSComputeRHSFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locG, void *user)
287: {
288: DM plex;
289: IS allcellIS;
290: PetscInt Nds, s;
292: PetscFunctionBegin;
293: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
294: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
295: PetscCall(DMGetNumDS(dm, &Nds));
296: for (s = 0; s < Nds; ++s) {
297: PetscDS ds;
298: IS cellIS;
299: PetscFormKey key;
301: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
302: key.value = 0;
303: key.field = 0;
304: key.part = 100;
305: if (!key.label) {
306: PetscCall(PetscObjectReference((PetscObject)allcellIS));
307: cellIS = allcellIS;
308: } else {
309: IS pointIS;
311: key.value = 1;
312: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
313: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
314: PetscCall(ISDestroy(&pointIS));
315: }
316: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locG, user));
317: PetscCall(ISDestroy(&cellIS));
318: }
319: PetscCall(ISDestroy(&allcellIS));
320: PetscCall(DMDestroy(&plex));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: /*@
325: DMTSCheckResidual - Check the residual of the exact solution
327: Input Parameters:
328: + ts - the `TS` object
329: . dm - the `DM`
330: . t - the time
331: . u - a `DM` vector
332: . u_t - a `DM` vector
333: - tol - A tolerance for the check, or -1 to print the results instead
335: Output Parameter:
336: . residual - The residual norm of the exact solution, or `NULL`
338: Level: developer
340: .seealso: [](ch_ts), `DM`, `DMTSCheckFromOptions()`, `DMTSCheckJacobian()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
341: @*/
342: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
343: {
344: MPI_Comm comm;
345: Vec r;
346: PetscReal res;
348: PetscFunctionBegin;
352: if (residual) PetscAssertPointer(residual, 7);
353: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
354: PetscCall(DMComputeExactSolution(dm, t, u, u_t));
355: PetscCall(VecDuplicate(u, &r));
356: PetscCall(TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE));
357: PetscCall(VecNorm(r, NORM_2, &res));
358: if (tol >= 0.0) {
359: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
360: } else if (residual) {
361: *residual = res;
362: } else {
363: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
364: PetscCall(VecFilter(r, 1.0e-10));
365: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)dm));
366: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
367: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
368: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
369: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL));
370: }
371: PetscCall(VecDestroy(&r));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /*@
376: DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
378: Input Parameters:
379: + ts - the `TS` object
380: . dm - the `DM`
381: . t - the time
382: . u - a `DM` vector
383: . u_t - a `DM` vector
384: - tol - A tolerance for the check, or -1 to print the results instead
386: Output Parameters:
387: + isLinear - Flag indicaing that the function looks linear, or `NULL`
388: - convRate - The rate of convergence of the linear model, or `NULL`
390: Level: developer
392: .seealso: [](ch_ts), `DNTSCheckFromOptions()`, `DMTSCheckResidual()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
393: @*/
394: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
395: {
396: MPI_Comm comm;
397: PetscDS ds;
398: Mat J, M;
399: MatNullSpace nullspace;
400: PetscReal dt, shift, slope, intercept;
401: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
403: PetscFunctionBegin;
407: if (isLinear) PetscAssertPointer(isLinear, 7);
408: if (convRate) PetscAssertPointer(convRate, 8);
409: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
410: PetscCall(DMComputeExactSolution(dm, t, u, u_t));
411: /* Create and view matrices */
412: PetscCall(TSGetTimeStep(ts, &dt));
413: shift = 1.0 / dt;
414: PetscCall(DMCreateMatrix(dm, &J));
415: PetscCall(DMGetDS(dm, &ds));
416: PetscCall(PetscDSHasJacobian(ds, &hasJac));
417: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
418: if (hasJac && hasPrec) {
419: PetscCall(DMCreateMatrix(dm, &M));
420: PetscCall(TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE));
421: PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
422: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
423: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
424: PetscCall(MatDestroy(&M));
425: } else {
426: PetscCall(TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE));
427: }
428: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
429: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
430: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
431: /* Check nullspace */
432: PetscCall(MatGetNullSpace(J, &nullspace));
433: if (nullspace) {
434: PetscBool isNull;
435: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
436: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
437: }
438: /* Taylor test */
439: {
440: PetscRandom rand;
441: Vec du, uhat, uhat_t, r, rhat, df;
442: PetscReal h;
443: PetscReal *es, *hs, *errors;
444: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
445: PetscInt Nv, v;
447: /* Choose a perturbation direction */
448: PetscCall(PetscRandomCreate(comm, &rand));
449: PetscCall(VecDuplicate(u, &du));
450: PetscCall(VecSetRandom(du, rand));
451: PetscCall(PetscRandomDestroy(&rand));
452: PetscCall(VecDuplicate(u, &df));
453: PetscCall(MatMult(J, du, df));
454: /* Evaluate residual at u, F(u), save in vector r */
455: PetscCall(VecDuplicate(u, &r));
456: PetscCall(TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE));
457: /* Look at the convergence of our Taylor approximation as we approach u */
458: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
459: PetscCall(PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors));
460: PetscCall(VecDuplicate(u, &uhat));
461: PetscCall(VecDuplicate(u, &uhat_t));
462: PetscCall(VecDuplicate(u, &rhat));
463: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
464: PetscCall(VecWAXPY(uhat, h, du, u));
465: PetscCall(VecWAXPY(uhat_t, h * shift, du, u_t));
466: /* 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 */
467: PetscCall(TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE));
468: PetscCall(VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df));
469: PetscCall(VecNorm(rhat, NORM_2, &errors[Nv]));
471: es[Nv] = PetscLog10Real(errors[Nv]);
472: hs[Nv] = PetscLog10Real(h);
473: }
474: PetscCall(VecDestroy(&uhat));
475: PetscCall(VecDestroy(&uhat_t));
476: PetscCall(VecDestroy(&rhat));
477: PetscCall(VecDestroy(&df));
478: PetscCall(VecDestroy(&r));
479: PetscCall(VecDestroy(&du));
480: for (v = 0; v < Nv; ++v) {
481: if ((tol >= 0) && (errors[v] > tol)) break;
482: else if (errors[v] > PETSC_SMALL) break;
483: }
484: if (v == Nv) isLin = PETSC_TRUE;
485: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
486: PetscCall(PetscFree3(es, hs, errors));
487: /* Slope should be about 2 */
488: if (tol >= 0) {
489: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
490: } else if (isLinear || convRate) {
491: if (isLinear) *isLinear = isLin;
492: if (convRate) *convRate = slope;
493: } else {
494: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
495: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
496: }
497: }
498: PetscCall(MatDestroy(&J));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /*@C
503: DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information based on
504: values in the options database
506: Input Parameters:
507: + ts - the `TS` object
508: - u - representative `TS` vector
510: Level: developer
512: Note:
513: The user must call `PetscDSSetExactSolution()` beforehand
515: Developer Notes:
516: What is the purpose of `u`, does it need to already have a solution or some other value in it?
518: .seealso: `DMTS`
519: @*/
520: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
521: {
522: DM dm;
523: SNES snes;
524: Vec sol, u_t;
525: PetscReal t;
526: PetscBool check;
528: PetscFunctionBegin;
529: PetscCall(PetscOptionsHasName(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-dmts_check", &check));
530: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
531: PetscCall(VecDuplicate(u, &sol));
532: PetscCall(VecCopy(u, sol));
533: PetscCall(TSSetSolution(ts, u));
534: PetscCall(TSGetDM(ts, &dm));
535: PetscCall(TSSetUp(ts));
536: PetscCall(TSGetSNES(ts, &snes));
537: PetscCall(SNESSetSolution(snes, u));
539: PetscCall(TSGetTime(ts, &t));
540: PetscCall(DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL));
541: PetscCall(DMGetGlobalVector(dm, &u_t));
542: PetscCall(DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL));
543: PetscCall(DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL));
544: PetscCall(DMRestoreGlobalVector(dm, &u_t));
546: PetscCall(VecDestroy(&sol));
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }