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: if (copy) {
22: PetscCall(DMCopyDMTS(dm, *plex));
23: PetscCall(DMCopyDMSNES(dm, *plex));
24: PetscCall(DMCopyAuxiliaryVec(dm, *plex));
25: }
26: } else {
27: PetscCall(PetscObjectReference((PetscObject)*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: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@
74: 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: Input Parameters:
77: + dm - The mesh
78: . time - The time
79: . locX - Local solution
80: . locX_t - Local solution time derivative, or `NULL`
81: - user - The user context
83: Level: developer
85: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
86: @*/
87: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
88: {
89: DM plex;
90: Vec faceGeometryFVM = NULL;
91: PetscInt Nf, f;
93: PetscFunctionBegin;
94: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
95: PetscCall(DMGetNumFields(plex, &Nf));
96: if (!locX_t) {
97: /* This is the RHS part */
98: for (f = 0; f < Nf; f++) {
99: PetscObject obj;
100: PetscClassId id;
102: PetscCall(DMGetField(plex, f, NULL, &obj));
103: PetscCall(PetscObjectGetClassId(obj, &id));
104: if (id == PETSCFV_CLASSID) {
105: PetscCall(DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL));
106: break;
107: }
108: }
109: }
110: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL));
111: PetscCall(DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL));
112: PetscCall(DMDestroy(&plex));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@
117: DMPlexTSComputeIFunctionFEM - Form the local residual `locF` from the local input `locX` using pointwise functions specified by the user
119: Input Parameters:
120: + dm - The mesh
121: . time - The time
122: . locX - Local solution
123: . locX_t - Local solution time derivative, or `NULL`
124: - user - The user context
126: Output Parameter:
127: . locF - Local output vector
129: Level: developer
131: .seealso: [](ch_ts), `DMPLEX`, `TS`, `DMPlexTSComputeRHSFunctionFEM()`
132: @*/
133: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
134: {
135: DM plex;
136: IS allcellIS;
137: PetscInt Nds, s;
139: PetscFunctionBegin;
140: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
141: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
142: PetscCall(DMGetNumDS(dm, &Nds));
143: for (s = 0; s < Nds; ++s) {
144: PetscDS ds;
145: IS cellIS;
146: PetscFormKey key;
148: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
149: key.value = 0;
150: key.field = 0;
151: key.part = 0;
152: if (!key.label) {
153: PetscCall(PetscObjectReference((PetscObject)allcellIS));
154: cellIS = allcellIS;
155: } else {
156: IS pointIS;
158: key.value = 1;
159: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
160: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
161: PetscCall(ISDestroy(&pointIS));
162: }
163: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user));
164: PetscCall(ISDestroy(&cellIS));
165: }
166: PetscCall(ISDestroy(&allcellIS));
167: PetscCall(DMDestroy(&plex));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: /*@
172: DMPlexTSComputeIJacobianFEM - Form the Jacobian `Jac` from the local input `locX` using pointwise functions specified by the user
174: Input Parameters:
175: + dm - The mesh
176: . time - The time
177: . locX - Local solution
178: . locX_t - Local solution time derivative, or `NULL`
179: . X_tShift - The multiplicative parameter for dF/du_t
180: - user - The user context
182: Output Parameters:
183: + Jac - the Jacobian
184: - JacP - an additional approximation for the Jacobian to be used to compute the preconditioner (often is `Jac`)
186: Level: developer
188: .seealso: [](ch_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
189: @*/
190: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
191: {
192: DM plex;
193: IS allcellIS;
194: PetscBool hasJac, hasPrec;
195: PetscInt Nds, s;
197: PetscFunctionBegin;
198: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
199: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
200: PetscCall(DMGetNumDS(dm, &Nds));
201: for (s = 0; s < Nds; ++s) {
202: PetscDS ds;
203: IS cellIS;
204: PetscFormKey key;
206: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
207: key.value = 0;
208: key.field = 0;
209: key.part = 0;
210: if (!key.label) {
211: PetscCall(PetscObjectReference((PetscObject)allcellIS));
212: cellIS = allcellIS;
213: } else {
214: IS pointIS;
216: key.value = 1;
217: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
218: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
219: PetscCall(ISDestroy(&pointIS));
220: }
221: if (!s) {
222: PetscCall(PetscDSHasJacobian(ds, &hasJac));
223: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
224: if (hasJac && hasPrec) PetscCall(MatZeroEntries(Jac));
225: PetscCall(MatZeroEntries(JacP));
226: }
227: PetscCall(DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user));
228: PetscCall(ISDestroy(&cellIS));
229: }
230: PetscCall(ISDestroy(&allcellIS));
231: PetscCall(DMDestroy(&plex));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@
236: DMPlexTSComputeRHSFunctionFEM - Form the local residual `locG` from the local input `locX` using pointwise functions specified by the user
238: Input Parameters:
239: + dm - The mesh
240: . time - The time
241: . locX - Local solution
242: - user - The user context
244: Output Parameter:
245: . locG - Local output vector
247: Level: developer
249: .seealso: [](ch_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeIJacobianFEM()`
250: @*/
251: PetscErrorCode DMPlexTSComputeRHSFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locG, void *user)
252: {
253: DM plex;
254: IS allcellIS;
255: PetscInt Nds, s;
257: PetscFunctionBegin;
258: PetscCall(DMTSConvertPlex(dm, &plex, PETSC_TRUE));
259: PetscCall(DMPlexGetAllCells_Internal(plex, &allcellIS));
260: PetscCall(DMGetNumDS(dm, &Nds));
261: for (s = 0; s < Nds; ++s) {
262: PetscDS ds;
263: IS cellIS;
264: PetscFormKey key;
266: PetscCall(DMGetRegionNumDS(dm, s, &key.label, NULL, &ds, NULL));
267: key.value = 0;
268: key.field = 0;
269: key.part = 100;
270: if (!key.label) {
271: PetscCall(PetscObjectReference((PetscObject)allcellIS));
272: cellIS = allcellIS;
273: } else {
274: IS pointIS;
276: key.value = 1;
277: PetscCall(DMLabelGetStratumIS(key.label, key.value, &pointIS));
278: PetscCall(ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS));
279: PetscCall(ISDestroy(&pointIS));
280: }
281: PetscCall(DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locG, user));
282: PetscCall(ISDestroy(&cellIS));
283: }
284: PetscCall(ISDestroy(&allcellIS));
285: PetscCall(DMDestroy(&plex));
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*@C
290: DMTSCheckResidual - Check the residual of the exact solution
292: Input Parameters:
293: + ts - the `TS` object
294: . dm - the `DM`
295: . t - the time
296: . u - a `DM` vector
297: . u_t - a `DM` vector
298: - tol - A tolerance for the check, or -1 to print the results instead
300: Output Parameter:
301: . residual - The residual norm of the exact solution, or `NULL`
303: Level: developer
305: .seealso: [](ch_ts), `DM`, `DMTSCheckFromOptions()`, `DMTSCheckJacobian()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
306: @*/
307: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
308: {
309: MPI_Comm comm;
310: Vec r;
311: PetscReal res;
313: PetscFunctionBegin;
317: if (residual) PetscAssertPointer(residual, 7);
318: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
319: PetscCall(DMComputeExactSolution(dm, t, u, u_t));
320: PetscCall(VecDuplicate(u, &r));
321: PetscCall(TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE));
322: PetscCall(VecNorm(r, NORM_2, &res));
323: if (tol >= 0.0) {
324: PetscCheck(res <= tol, comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double)res, (double)tol);
325: } else if (residual) {
326: *residual = res;
327: } else {
328: PetscCall(PetscPrintf(comm, "L_2 Residual: %g\n", (double)res));
329: PetscCall(VecFilter(r, 1.0e-10));
330: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)dm));
331: PetscCall(PetscObjectSetName((PetscObject)r, "Initial Residual"));
332: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)r, "res_"));
333: PetscCall(VecViewFromOptions(r, NULL, "-vec_view"));
334: PetscCall(PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL));
335: }
336: PetscCall(VecDestroy(&r));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@C
341: DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
343: Input Parameters:
344: + ts - the `TS` object
345: . dm - the `DM`
346: . t - the time
347: . u - a `DM` vector
348: . u_t - a `DM` vector
349: - tol - A tolerance for the check, or -1 to print the results instead
351: Output Parameters:
352: + isLinear - Flag indicaing that the function looks linear, or `NULL`
353: - convRate - The rate of convergence of the linear model, or `NULL`
355: Level: developer
357: .seealso: [](ch_ts), `DNTSCheckFromOptions()`, `DMTSCheckResidual()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
358: @*/
359: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
360: {
361: MPI_Comm comm;
362: PetscDS ds;
363: Mat J, M;
364: MatNullSpace nullspace;
365: PetscReal dt, shift, slope, intercept;
366: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
368: PetscFunctionBegin;
372: if (isLinear) PetscAssertPointer(isLinear, 7);
373: if (convRate) PetscAssertPointer(convRate, 8);
374: PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
375: PetscCall(DMComputeExactSolution(dm, t, u, u_t));
376: /* Create and view matrices */
377: PetscCall(TSGetTimeStep(ts, &dt));
378: shift = 1.0 / dt;
379: PetscCall(DMCreateMatrix(dm, &J));
380: PetscCall(DMGetDS(dm, &ds));
381: PetscCall(PetscDSHasJacobian(ds, &hasJac));
382: PetscCall(PetscDSHasJacobianPreconditioner(ds, &hasPrec));
383: if (hasJac && hasPrec) {
384: PetscCall(DMCreateMatrix(dm, &M));
385: PetscCall(TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE));
386: PetscCall(PetscObjectSetName((PetscObject)M, "Preconditioning Matrix"));
387: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_"));
388: PetscCall(MatViewFromOptions(M, NULL, "-mat_view"));
389: PetscCall(MatDestroy(&M));
390: } else {
391: PetscCall(TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE));
392: }
393: PetscCall(PetscObjectSetName((PetscObject)J, "Jacobian"));
394: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)J, "jac_"));
395: PetscCall(MatViewFromOptions(J, NULL, "-mat_view"));
396: /* Check nullspace */
397: PetscCall(MatGetNullSpace(J, &nullspace));
398: if (nullspace) {
399: PetscBool isNull;
400: PetscCall(MatNullSpaceTest(nullspace, J, &isNull));
401: PetscCheck(isNull, comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
402: }
403: /* Taylor test */
404: {
405: PetscRandom rand;
406: Vec du, uhat, uhat_t, r, rhat, df;
407: PetscReal h;
408: PetscReal *es, *hs, *errors;
409: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
410: PetscInt Nv, v;
412: /* Choose a perturbation direction */
413: PetscCall(PetscRandomCreate(comm, &rand));
414: PetscCall(VecDuplicate(u, &du));
415: PetscCall(VecSetRandom(du, rand));
416: PetscCall(PetscRandomDestroy(&rand));
417: PetscCall(VecDuplicate(u, &df));
418: PetscCall(MatMult(J, du, df));
419: /* Evaluate residual at u, F(u), save in vector r */
420: PetscCall(VecDuplicate(u, &r));
421: PetscCall(TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE));
422: /* Look at the convergence of our Taylor approximation as we approach u */
423: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
424: ;
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: }