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