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: /*@
34: DMPlexTSComputeRHSFunctionFVM - Form the forcing `F` from the local input `locX` using pointwise functions specified by the user
36: Input Parameters:
37: + dm - The mesh
38: . time - The time
39: . locX - Local solution
40: - user - The user context
42: Output Parameter:
43: . F - Global output vector
45: Level: developer
47: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
48: @*/
49: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
50: {
51: Vec locF;
52: IS cellIS;
53: DM plex;
54: PetscInt depth;
55: PetscFormKey key = {NULL, 0, 0, 0};
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(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locF, user));
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: }
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
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: - user - The user context
84: Level: developer
86: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
87: @*/
88: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
89: {
90: DM plex;
91: Vec faceGeometryFVM = NULL;
92: PetscInt Nf, f;
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;
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: }
117: /*@
118: DMPlexTSComputeIFunctionFEM - Form the local residual `locF` from the local input `locX` using pointwise functions specified by the user
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: - user - The user context
127: Output Parameter:
128: . locF - Local output vector
130: Level: developer
132: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexTSComputeRHSFunctionFEM()`
133: @*/
134: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
135: {
136: DM plex;
137: IS allcellIS;
138: PetscInt Nds, s;
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;
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;
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(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user));
165: PetscCall(ISDestroy(&cellIS));
166: }
167: PetscCall(ISDestroy(&allcellIS));
168: PetscCall(DMDestroy(&plex));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: /*@
173: DMPlexTSComputeIJacobianFEM - Form the Jacobian `Jac` from the local input `locX` using pointwise functions specified by the user
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: - user - The user context
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`)
187: Level: developer
189: .seealso: [](ch_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
190: @*/
191: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
192: {
193: DM plex;
194: IS allcellIS;
195: PetscBool hasJac, hasPrec;
196: PetscInt Nds, s;
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;
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;
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(DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user));
229: PetscCall(ISDestroy(&cellIS));
230: }
231: PetscCall(ISDestroy(&allcellIS));
232: PetscCall(DMDestroy(&plex));
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
236: /*@
237: DMPlexTSComputeRHSFunctionFEM - Form the local residual `locG` from the local input `locX` using pointwise functions specified by the user
239: Input Parameters:
240: + dm - The mesh
241: . time - The time
242: . locX - Local solution
243: - user - The user context
245: Output Parameter:
246: . locG - Local output vector
248: Level: developer
250: .seealso: [](ch_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeIJacobianFEM()`
251: @*/
252: PetscErrorCode DMPlexTSComputeRHSFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locG, void *user)
253: {
254: DM plex;
255: IS allcellIS;
256: PetscInt Nds, s;
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;
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;
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(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locG, user));
283: PetscCall(ISDestroy(&cellIS));
284: }
285: PetscCall(ISDestroy(&allcellIS));
286: PetscCall(DMDestroy(&plex));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /*@
291: DMTSCheckResidual - Check the residual of the exact solution
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
301: Output Parameter:
302: . residual - The residual norm of the exact solution, or `NULL`
304: Level: developer
306: .seealso: [](ch_ts), `DM`, `DMTSCheckFromOptions()`, `DMTSCheckJacobian()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
307: @*/
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;
314: PetscFunctionBegin;
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: }
341: /*@
342: DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
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
352: Output Parameters:
353: + isLinear - Flag indicaing that the function looks linear, or `NULL`
354: - convRate - The rate of convergence of the linear model, or `NULL`
356: Level: developer
358: .seealso: [](ch_ts), `DNTSCheckFromOptions()`, `DMTSCheckResidual()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
359: @*/
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;
369: PetscFunctionBegin;
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, "Preconditioning Matrix"));
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;
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]));
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) && (errors[v] > tol)) break;
448: else if (errors[v] > PETSC_SMALL) break;
449: }
450: if (v == Nv) isLin = PETSC_TRUE;
451: PetscCall(PetscLinearRegression(Nv, hs, es, &slope, &intercept));
452: PetscCall(PetscFree3(es, hs, errors));
453: /* Slope should be about 2 */
454: if (tol >= 0) {
455: PetscCheck(isLin || PetscAbsReal(2 - slope) <= tol, comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double)slope);
456: } else if (isLinear || convRate) {
457: if (isLinear) *isLinear = isLin;
458: if (convRate) *convRate = slope;
459: } else {
460: if (!isLin) PetscCall(PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope));
461: else PetscCall(PetscPrintf(comm, "Function appears to be linear\n"));
462: }
463: }
464: PetscCall(MatDestroy(&J));
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@C
469: DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information based on
470: values in the options database
472: Input Parameters:
473: + ts - the `TS` object
474: - u - representative `TS` vector
476: Level: developer
478: Note:
479: The user must call `PetscDSSetExactSolution()` beforehand
481: Developer Notes:
482: What is the purpose of `u`, does it need to already have a solution or some other value in it?
484: .seealso: `DMTS`
485: @*/
486: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
487: {
488: DM dm;
489: SNES snes;
490: Vec sol, u_t;
491: PetscReal t;
492: PetscBool check;
494: PetscFunctionBegin;
495: PetscCall(PetscOptionsHasName(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-dmts_check", &check));
496: if (!check) PetscFunctionReturn(PETSC_SUCCESS);
497: PetscCall(VecDuplicate(u, &sol));
498: PetscCall(VecCopy(u, sol));
499: PetscCall(TSSetSolution(ts, u));
500: PetscCall(TSGetDM(ts, &dm));
501: PetscCall(TSSetUp(ts));
502: PetscCall(TSGetSNES(ts, &snes));
503: PetscCall(SNESSetSolution(snes, u));
505: PetscCall(TSGetTime(ts, &t));
506: PetscCall(DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL));
507: PetscCall(DMGetGlobalVector(dm, &u_t));
508: PetscCall(DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL));
509: PetscCall(DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL));
510: PetscCall(DMRestoreGlobalVector(dm, &u_t));
512: PetscCall(VecDestroy(&sol));
513: PetscFunctionReturn(PETSC_SUCCESS);
514: }