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: PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isPlex);
12: if (isPlex) {
13: *plex = dm;
14: PetscObjectReference((PetscObject)dm);
15: } else {
16: PetscObjectQuery((PetscObject)dm, "dm_plex", (PetscObject *)plex);
17: if (!*plex) {
18: DMConvert(dm, DMPLEX, plex);
19: PetscObjectCompose((PetscObject)dm, "dm_plex", (PetscObject)*plex);
20: if (copy) {
21: DMCopyDMTS(dm, *plex);
22: DMCopyDMSNES(dm, *plex);
23: DMCopyAuxiliaryVec(dm, *plex);
24: }
25: } else {
26: PetscObjectReference((PetscObject)*plex);
27: }
28: }
29: return 0;
30: }
32: /*@
33: DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user
35: Input Parameters:
36: + dm - The mesh
37: . t - The time
38: . locX - Local solution
39: - user - The user context
41: Output Parameter:
42: . F - Global output vector
44: Level: developer
46: .seealso: [](chapter_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
47: @*/
48: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
49: {
50: Vec locF;
51: IS cellIS;
52: DM plex;
53: PetscInt depth;
54: PetscFormKey key = {NULL, 0, 0, 0};
56: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
57: DMPlexGetDepth(plex, &depth);
58: DMGetStratumIS(plex, "dim", depth, &cellIS);
59: if (!cellIS) DMGetStratumIS(plex, "depth", depth, &cellIS);
60: DMGetLocalVector(plex, &locF);
61: VecZeroEntries(locF);
62: DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locF, user);
63: DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
64: DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
65: DMRestoreLocalVector(plex, &locF);
66: ISDestroy(&cellIS);
67: DMDestroy(&plex);
68: return 0;
69: }
71: /*@
72: 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
74: Input Parameters:
75: + dm - The mesh
76: . t - The time
77: . locX - Local solution
78: . locX_t - Local solution time derivative, or NULL
79: - user - The user context
81: Level: developer
83: .seealso: [](chapter_ts), `DMPLEX`, `TS`, `DMPlexComputeJacobianActionFEM()`
84: @*/
85: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
86: {
87: DM plex;
88: Vec faceGeometryFVM = NULL;
89: PetscInt Nf, f;
91: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
92: DMGetNumFields(plex, &Nf);
93: if (!locX_t) {
94: /* This is the RHS part */
95: for (f = 0; f < Nf; f++) {
96: PetscObject obj;
97: PetscClassId id;
99: DMGetField(plex, f, NULL, &obj);
100: PetscObjectGetClassId(obj, &id);
101: if (id == PETSCFV_CLASSID) {
102: DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
103: break;
104: }
105: }
106: }
107: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
108: DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
109: DMDestroy(&plex);
110: return 0;
111: }
113: /*@
114: DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
116: Input Parameters:
117: + dm - The mesh
118: . t - The time
119: . locX - Local solution
120: . locX_t - Local solution time derivative, or NULL
121: - user - The user context
123: Output Parameter:
124: . locF - Local output vector
126: Level: developer
128: .seealso: [](chapter_ts), `DMPLEX`, `TS`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
129: @*/
130: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
131: {
132: DM plex;
133: IS allcellIS;
134: PetscInt Nds, s;
136: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
137: DMPlexGetAllCells_Internal(plex, &allcellIS);
138: DMGetNumDS(dm, &Nds);
139: for (s = 0; s < Nds; ++s) {
140: PetscDS ds;
141: IS cellIS;
142: PetscFormKey key;
144: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
145: key.value = 0;
146: key.field = 0;
147: key.part = 0;
148: if (!key.label) {
149: PetscObjectReference((PetscObject)allcellIS);
150: cellIS = allcellIS;
151: } else {
152: IS pointIS;
154: key.value = 1;
155: DMLabelGetStratumIS(key.label, key.value, &pointIS);
156: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
157: ISDestroy(&pointIS);
158: }
159: DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user);
160: ISDestroy(&cellIS);
161: }
162: ISDestroy(&allcellIS);
163: DMDestroy(&plex);
164: return 0;
165: }
167: /*@
168: DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
170: Input Parameters:
171: + dm - The mesh
172: . t - The time
173: . locX - Local solution
174: . locX_t - Local solution time derivative, or NULL
175: . X_tshift - The multiplicative parameter for dF/du_t
176: - user - The user context
178: Output Parameter:
179: . locF - Local output vector
181: Level: developer
183: .seealso: [](chapter_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeRHSFunctionFEM()`
184: @*/
185: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
186: {
187: DM plex;
188: IS allcellIS;
189: PetscBool hasJac, hasPrec;
190: PetscInt Nds, s;
192: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
193: DMPlexGetAllCells_Internal(plex, &allcellIS);
194: DMGetNumDS(dm, &Nds);
195: for (s = 0; s < Nds; ++s) {
196: PetscDS ds;
197: IS cellIS;
198: PetscFormKey key;
200: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
201: key.value = 0;
202: key.field = 0;
203: key.part = 0;
204: if (!key.label) {
205: PetscObjectReference((PetscObject)allcellIS);
206: cellIS = allcellIS;
207: } else {
208: IS pointIS;
210: key.value = 1;
211: DMLabelGetStratumIS(key.label, key.value, &pointIS);
212: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
213: ISDestroy(&pointIS);
214: }
215: if (!s) {
216: PetscDSHasJacobian(ds, &hasJac);
217: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
218: if (hasJac && hasPrec) MatZeroEntries(Jac);
219: MatZeroEntries(JacP);
220: }
221: DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
222: ISDestroy(&cellIS);
223: }
224: ISDestroy(&allcellIS);
225: DMDestroy(&plex);
226: return 0;
227: }
229: /*@
230: DMPlexTSComputeRHSFunctionFEM - Form the local residual G from the local input X using pointwise functions specified by the user
232: Input Parameters:
233: + dm - The mesh
234: . t - The time
235: . locX - Local solution
236: - user - The user context
238: Output Parameter:
239: . locG - Local output vector
241: Level: developer
243: .seealso: [](chapter_ts), `TS`, `DM`, `DMPlexTSComputeIFunctionFEM()`, `DMPlexTSComputeIJacobianFEM()`
244: @*/
245: PetscErrorCode DMPlexTSComputeRHSFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locG, void *user)
246: {
247: DM plex;
248: IS allcellIS;
249: PetscInt Nds, s;
251: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
252: DMPlexGetAllCells_Internal(plex, &allcellIS);
253: DMGetNumDS(dm, &Nds);
254: for (s = 0; s < Nds; ++s) {
255: PetscDS ds;
256: IS cellIS;
257: PetscFormKey key;
259: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
260: key.value = 0;
261: key.field = 0;
262: key.part = 100;
263: if (!key.label) {
264: PetscObjectReference((PetscObject)allcellIS);
265: cellIS = allcellIS;
266: } else {
267: IS pointIS;
269: key.value = 1;
270: DMLabelGetStratumIS(key.label, key.value, &pointIS);
271: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
272: ISDestroy(&pointIS);
273: }
274: DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locG, user);
275: ISDestroy(&cellIS);
276: }
277: ISDestroy(&allcellIS);
278: DMDestroy(&plex);
279: return 0;
280: }
282: /*@C
283: DMTSCheckResidual - Check the residual of the exact solution
285: Input Parameters:
286: + ts - the `TS` object
287: . dm - the `DM`
288: . t - the time
289: . u - a `DM` vector
290: . u_t - a `DM` vector
291: - tol - A tolerance for the check, or -1 to print the results instead
293: Output Parameters:
294: . residual - The residual norm of the exact solution, or NULL
296: Level: developer
298: .seealso: [](chapter_ts), `DM`, `DMTSCheckFromOptions()`, `DMTSCheckJacobian()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckJacobian()`
299: @*/
300: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
301: {
302: MPI_Comm comm;
303: Vec r;
304: PetscReal res;
310: PetscObjectGetComm((PetscObject)ts, &comm);
311: DMComputeExactSolution(dm, t, u, u_t);
312: VecDuplicate(u, &r);
313: TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
314: VecNorm(r, NORM_2, &res);
315: if (tol >= 0.0) {
317: } else if (residual) {
318: *residual = res;
319: } else {
320: PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
321: VecChop(r, 1.0e-10);
322: PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", (PetscObject)dm);
323: PetscObjectSetName((PetscObject)r, "Initial Residual");
324: PetscObjectSetOptionsPrefix((PetscObject)r, "res_");
325: VecViewFromOptions(r, NULL, "-vec_view");
326: PetscObjectCompose((PetscObject)r, "__Vec_bc_zero__", NULL);
327: }
328: VecDestroy(&r);
329: return 0;
330: }
332: /*@C
333: DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
335: Input Parameters:
336: + ts - the TS object
337: . dm - the DM
338: . t - the time
339: . u - a DM vector
340: . u_t - a DM vector
341: - tol - A tolerance for the check, or -1 to print the results instead
343: Output Parameters:
344: + isLinear - Flag indicaing that the function looks linear, or NULL
345: - convRate - The rate of convergence of the linear model, or NULL
347: Level: developer
349: .seealso: [](chapter_ts), `DNTSCheckFromOptions()`, `DMTSCheckResidual()`, `DNSNESCheckFromOptions()`, `DMSNESCheckDiscretization()`, `DMSNESCheckResidual()`
350: @*/
351: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
352: {
353: MPI_Comm comm;
354: PetscDS ds;
355: Mat J, M;
356: MatNullSpace nullspace;
357: PetscReal dt, shift, slope, intercept;
358: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
365: PetscObjectGetComm((PetscObject)ts, &comm);
366: DMComputeExactSolution(dm, t, u, u_t);
367: /* Create and view matrices */
368: TSGetTimeStep(ts, &dt);
369: shift = 1.0 / dt;
370: DMCreateMatrix(dm, &J);
371: DMGetDS(dm, &ds);
372: PetscDSHasJacobian(ds, &hasJac);
373: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
374: if (hasJac && hasPrec) {
375: DMCreateMatrix(dm, &M);
376: TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
377: PetscObjectSetName((PetscObject)M, "Preconditioning Matrix");
378: PetscObjectSetOptionsPrefix((PetscObject)M, "jacpre_");
379: MatViewFromOptions(M, NULL, "-mat_view");
380: MatDestroy(&M);
381: } else {
382: TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
383: }
384: PetscObjectSetName((PetscObject)J, "Jacobian");
385: PetscObjectSetOptionsPrefix((PetscObject)J, "jac_");
386: MatViewFromOptions(J, NULL, "-mat_view");
387: /* Check nullspace */
388: MatGetNullSpace(J, &nullspace);
389: if (nullspace) {
390: PetscBool isNull;
391: MatNullSpaceTest(nullspace, J, &isNull);
393: }
394: /* Taylor test */
395: {
396: PetscRandom rand;
397: Vec du, uhat, uhat_t, r, rhat, df;
398: PetscReal h;
399: PetscReal *es, *hs, *errors;
400: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
401: PetscInt Nv, v;
403: /* Choose a perturbation direction */
404: PetscRandomCreate(comm, &rand);
405: VecDuplicate(u, &du);
406: VecSetRandom(du, rand);
407: PetscRandomDestroy(&rand);
408: VecDuplicate(u, &df);
409: MatMult(J, du, df);
410: /* Evaluate residual at u, F(u), save in vector r */
411: VecDuplicate(u, &r);
412: TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
413: /* Look at the convergence of our Taylor approximation as we approach u */
414: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv)
415: ;
416: PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
417: VecDuplicate(u, &uhat);
418: VecDuplicate(u, &uhat_t);
419: VecDuplicate(u, &rhat);
420: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
421: VecWAXPY(uhat, h, du, u);
422: VecWAXPY(uhat_t, h * shift, du, u_t);
423: /* 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 */
424: TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
425: VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
426: VecNorm(rhat, NORM_2, &errors[Nv]);
428: es[Nv] = PetscLog10Real(errors[Nv]);
429: hs[Nv] = PetscLog10Real(h);
430: }
431: VecDestroy(&uhat);
432: VecDestroy(&uhat_t);
433: VecDestroy(&rhat);
434: VecDestroy(&df);
435: VecDestroy(&r);
436: VecDestroy(&du);
437: for (v = 0; v < Nv; ++v) {
438: if ((tol >= 0) && (errors[v] > tol)) break;
439: else if (errors[v] > PETSC_SMALL) break;
440: }
441: if (v == Nv) isLin = PETSC_TRUE;
442: PetscLinearRegression(Nv, hs, es, &slope, &intercept);
443: PetscFree3(es, hs, errors);
444: /* Slope should be about 2 */
445: if (tol >= 0) {
447: } else if (isLinear || convRate) {
448: if (isLinear) *isLinear = isLin;
449: if (convRate) *convRate = slope;
450: } else {
451: if (!isLin) PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double)slope);
452: else PetscPrintf(comm, "Function appears to be linear\n");
453: }
454: }
455: MatDestroy(&J);
456: return 0;
457: }
459: /*@C
460: DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
462: Input Parameters:
463: + ts - the `TS` object
464: - u - representative `TS` vector
466: Note:
467: The user must call `PetscDSSetExactSolution()` beforehand
469: Level: developer
470: @*/
471: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
472: {
473: DM dm;
474: SNES snes;
475: Vec sol, u_t;
476: PetscReal t;
477: PetscBool check;
479: PetscOptionsHasName(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-dmts_check", &check);
480: if (!check) return 0;
481: VecDuplicate(u, &sol);
482: VecCopy(u, sol);
483: TSSetSolution(ts, u);
484: TSGetDM(ts, &dm);
485: TSSetUp(ts);
486: TSGetSNES(ts, &snes);
487: SNESSetSolution(snes, u);
489: TSGetTime(ts, &t);
490: DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
491: DMGetGlobalVector(dm, &u_t);
492: DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
493: DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
494: DMRestoreGlobalVector(dm, &u_t);
496: VecDestroy(&sol);
497: return 0;
498: }