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