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;
13: PetscObjectTypeCompare((PetscObject) dm, DMPLEX, &isPlex);
14: if (isPlex) {
15: *plex = dm;
16: PetscObjectReference((PetscObject) dm);
17: } else {
18: PetscObjectQuery((PetscObject) dm, "dm_plex", (PetscObject *) plex);
19: if (!*plex) {
20: DMConvert(dm,DMPLEX,plex);
21: PetscObjectCompose((PetscObject) dm, "dm_plex", (PetscObject) *plex);
22: if (copy) {
23: PetscInt i;
24: PetscObject obj;
25: const char *comps[3] = {"A","dmAux","dmCh"};
27: DMCopyDMTS(dm, *plex);
28: DMCopyDMSNES(dm, *plex);
29: for (i = 0; i < 3; i++) {
30: PetscObjectQuery((PetscObject) dm, comps[i], &obj);
31: PetscObjectCompose((PetscObject) *plex, comps[i], obj);
32: }
33: }
34: } else {
35: PetscObjectReference((PetscObject) *plex);
36: }
37: }
38: return(0);
39: }
41: /*@
42: DMPlexTSComputeRHSFunctionFVM - Form the local forcing F from the local input X using pointwise functions specified by the user
44: Input Parameters:
45: + dm - The mesh
46: . t - The time
47: . locX - Local solution
48: - user - The user context
50: Output Parameter:
51: . F - Global output vector
53: Level: developer
55: .seealso: DMPlexComputeJacobianActionFEM()
56: @*/
57: PetscErrorCode DMPlexTSComputeRHSFunctionFVM(DM dm, PetscReal time, Vec locX, Vec F, void *user)
58: {
59: Vec locF;
60: IS cellIS;
61: DM plex;
62: PetscInt depth;
63: PetscHashFormKey key = {NULL, 0, 0};
67: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
68: DMPlexGetDepth(plex, &depth);
69: DMGetStratumIS(plex, "dim", depth, &cellIS);
70: if (!cellIS) {
71: DMGetStratumIS(plex, "depth", depth, &cellIS);
72: }
73: DMGetLocalVector(plex, &locF);
74: VecZeroEntries(locF);
75: DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, NULL, time, locF, user);
76: DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
77: DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
78: DMRestoreLocalVector(plex, &locF);
79: ISDestroy(&cellIS);
80: DMDestroy(&plex);
81: return(0);
82: }
84: /*@
85: 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
87: Input Parameters:
88: + dm - The mesh
89: . t - The time
90: . locX - Local solution
91: . locX_t - Local solution time derivative, or NULL
92: - user - The user context
94: Level: developer
96: .seealso: DMPlexComputeJacobianActionFEM()
97: @*/
98: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
99: {
100: DM plex;
101: Vec faceGeometryFVM = NULL;
102: PetscInt Nf, f;
106: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
107: DMGetNumFields(plex, &Nf);
108: if (!locX_t) {
109: /* This is the RHS part */
110: for (f = 0; f < Nf; f++) {
111: PetscObject obj;
112: PetscClassId id;
114: DMGetField(plex, f, NULL, &obj);
115: PetscObjectGetClassId(obj, &id);
116: if (id == PETSCFV_CLASSID) {
117: DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
118: break;
119: }
120: }
121: }
122: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
123: DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
124: DMDestroy(&plex);
125: return(0);
126: }
128: /*@
129: DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
131: Input Parameters:
132: + dm - The mesh
133: . t - The time
134: . locX - Local solution
135: . locX_t - Local solution time derivative, or NULL
136: - user - The user context
138: Output Parameter:
139: . locF - Local output vector
141: Level: developer
143: .seealso: DMPlexComputeJacobianActionFEM()
144: @*/
145: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
146: {
147: DM plex;
148: IS allcellIS;
149: PetscInt Nds, s;
153: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
154: DMPlexGetAllCells_Internal(plex, &allcellIS);
155: DMGetNumDS(dm, &Nds);
156: for (s = 0; s < Nds; ++s) {
157: PetscDS ds;
158: IS cellIS;
159: PetscHashFormKey key;
161: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
162: key.value = 0;
163: key.field = 0;
164: if (!key.label) {
165: PetscObjectReference((PetscObject) allcellIS);
166: cellIS = allcellIS;
167: } else {
168: IS pointIS;
170: key.value = 1;
171: DMLabelGetStratumIS(key.label, key.value, &pointIS);
172: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
173: ISDestroy(&pointIS);
174: }
175: DMPlexComputeResidual_Internal(plex, key, cellIS, time, locX, locX_t, time, locF, user);
176: ISDestroy(&cellIS);
177: }
178: ISDestroy(&allcellIS);
179: DMDestroy(&plex);
180: return(0);
181: }
183: /*@
184: DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
186: Input Parameters:
187: + dm - The mesh
188: . t - The time
189: . locX - Local solution
190: . locX_t - Local solution time derivative, or NULL
191: . X_tshift - The multiplicative parameter for dF/du_t
192: - user - The user context
194: Output Parameter:
195: . locF - Local output vector
197: Level: developer
199: .seealso: DMPlexComputeJacobianActionFEM()
200: @*/
201: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
202: {
203: DM plex;
204: IS allcellIS;
205: PetscBool hasJac, hasPrec;
206: PetscInt Nds, s;
210: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
211: DMPlexGetAllCells_Internal(plex, &allcellIS);
212: DMGetNumDS(dm, &Nds);
213: for (s = 0; s < Nds; ++s) {
214: PetscDS ds;
215: IS cellIS;
216: PetscHashFormKey key;
218: DMGetRegionNumDS(dm, s, &key.label, NULL, &ds);
219: key.value = 0;
220: key.field = 0;
221: if (!key.label) {
222: PetscObjectReference((PetscObject) allcellIS);
223: cellIS = allcellIS;
224: } else {
225: IS pointIS;
227: key.value = 1;
228: DMLabelGetStratumIS(key.label, key.value, &pointIS);
229: ISIntersect_Caching_Internal(allcellIS, pointIS, &cellIS);
230: ISDestroy(&pointIS);
231: }
232: if (!s) {
233: PetscDSHasJacobian(ds, &hasJac);
234: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
235: if (hasJac && hasPrec) {MatZeroEntries(Jac);}
236: MatZeroEntries(JacP);
237: }
238: DMPlexComputeJacobian_Internal(plex, key, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
239: ISDestroy(&cellIS);
240: }
241: ISDestroy(&allcellIS);
242: DMDestroy(&plex);
243: return(0);
244: }
246: /*@C
247: DMTSCheckResidual - Check the residual of the exact solution
249: Input Parameters:
250: + ts - the TS object
251: . dm - the DM
252: . t - the time
253: . u - a DM vector
254: . u_t - a DM vector
255: - tol - A tolerance for the check, or -1 to print the results instead
257: Output Parameters:
258: . residual - The residual norm of the exact solution, or NULL
260: Level: developer
262: .seealso: DNTSCheckFromOptions(), DMTSCheckJacobian(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
263: @*/
264: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
265: {
266: MPI_Comm comm;
267: Vec r;
268: PetscReal res;
276: PetscObjectGetComm((PetscObject) ts, &comm);
277: DMComputeExactSolution(dm, t, u, u_t);
278: VecDuplicate(u, &r);
279: TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
280: VecNorm(r, NORM_2, &res);
281: if (tol >= 0.0) {
282: if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
283: } else if (residual) {
284: *residual = res;
285: } else {
286: PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
287: VecChop(r, 1.0e-10);
288: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) dm);
289: PetscObjectSetName((PetscObject) r, "Initial Residual");
290: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
291: VecViewFromOptions(r, NULL, "-vec_view");
292: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
293: }
294: VecDestroy(&r);
295: return(0);
296: }
298: /*@C
299: DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
301: Input Parameters:
302: + ts - the TS object
303: . dm - the DM
304: . t - the time
305: . u - a DM vector
306: . u_t - a DM vector
307: - tol - A tolerance for the check, or -1 to print the results instead
309: Output Parameters:
310: + isLinear - Flag indicaing that the function looks linear, or NULL
311: - convRate - The rate of convergence of the linear model, or NULL
313: Level: developer
315: .seealso: DNTSCheckFromOptions(), DMTSCheckResidual(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
316: @*/
317: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
318: {
319: MPI_Comm comm;
320: PetscDS ds;
321: Mat J, M;
322: MatNullSpace nullspace;
323: PetscReal dt, shift, slope, intercept;
324: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
333: PetscObjectGetComm((PetscObject) ts, &comm);
334: DMComputeExactSolution(dm, t, u, u_t);
335: /* Create and view matrices */
336: TSGetTimeStep(ts, &dt);
337: shift = 1.0/dt;
338: DMCreateMatrix(dm, &J);
339: DMGetDS(dm, &ds);
340: PetscDSHasJacobian(ds, &hasJac);
341: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
342: if (hasJac && hasPrec) {
343: DMCreateMatrix(dm, &M);
344: TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
345: PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
346: PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
347: MatViewFromOptions(M, NULL, "-mat_view");
348: MatDestroy(&M);
349: } else {
350: TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
351: }
352: PetscObjectSetName((PetscObject) J, "Jacobian");
353: PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
354: MatViewFromOptions(J, NULL, "-mat_view");
355: /* Check nullspace */
356: MatGetNullSpace(J, &nullspace);
357: if (nullspace) {
358: PetscBool isNull;
359: MatNullSpaceTest(nullspace, J, &isNull);
360: if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
361: }
362: /* Taylor test */
363: {
364: PetscRandom rand;
365: Vec du, uhat, uhat_t, r, rhat, df;
366: PetscReal h;
367: PetscReal *es, *hs, *errors;
368: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
369: PetscInt Nv, v;
371: /* Choose a perturbation direction */
372: PetscRandomCreate(comm, &rand);
373: VecDuplicate(u, &du);
374: VecSetRandom(du, rand);
375: PetscRandomDestroy(&rand);
376: VecDuplicate(u, &df);
377: MatMult(J, du, df);
378: /* Evaluate residual at u, F(u), save in vector r */
379: VecDuplicate(u, &r);
380: TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
381: /* Look at the convergence of our Taylor approximation as we approach u */
382: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
383: PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
384: VecDuplicate(u, &uhat);
385: VecDuplicate(u, &uhat_t);
386: VecDuplicate(u, &rhat);
387: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
388: VecWAXPY(uhat, h, du, u);
389: VecWAXPY(uhat_t, h*shift, du, u_t);
390: /* 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 */
391: TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
392: VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
393: VecNorm(rhat, NORM_2, &errors[Nv]);
395: es[Nv] = PetscLog10Real(errors[Nv]);
396: hs[Nv] = PetscLog10Real(h);
397: }
398: VecDestroy(&uhat);
399: VecDestroy(&uhat_t);
400: VecDestroy(&rhat);
401: VecDestroy(&df);
402: VecDestroy(&r);
403: VecDestroy(&du);
404: for (v = 0; v < Nv; ++v) {
405: if ((tol >= 0) && (errors[v] > tol)) break;
406: else if (errors[v] > PETSC_SMALL) break;
407: }
408: if (v == Nv) isLin = PETSC_TRUE;
409: PetscLinearRegression(Nv, hs, es, &slope, &intercept);
410: PetscFree3(es, hs, errors);
411: /* Slope should be about 2 */
412: if (tol >= 0) {
413: if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
414: } else if (isLinear || convRate) {
415: if (isLinear) *isLinear = isLin;
416: if (convRate) *convRate = slope;
417: } else {
418: if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
419: else {PetscPrintf(comm, "Function appears to be linear\n");}
420: }
421: }
422: MatDestroy(&J);
423: return(0);
424: }
426: /*@C
427: DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
429: Input Parameters:
430: + ts - the TS object
431: - u - representative TS vector
433: Note: The user must call PetscDSSetExactSolution() beforehand
435: Level: developer
436: @*/
437: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
438: {
439: DM dm;
440: SNES snes;
441: Vec sol, u_t;
442: PetscReal t;
443: PetscBool check;
447: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
448: if (!check) return(0);
449: VecDuplicate(u, &sol);
450: VecCopy(u, sol);
451: TSSetSolution(ts, u);
452: TSGetDM(ts, &dm);
453: TSSetUp(ts);
454: TSGetSNES(ts, &snes);
455: SNESSetSolution(snes, u);
457: TSGetTime(ts, &t);
458: DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
459: DMGetGlobalVector(dm, &u_t);
460: DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
461: DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
462: DMRestoreGlobalVector(dm, &u_t);
464: VecDestroy(&sol);
465: return(0);
466: }