Actual source code: dmplexts.c
petsc-3.14.6 2021-03-30
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;
66: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
67: DMPlexGetDepth(plex, &depth);
68: DMGetStratumIS(plex, "dim", depth, &cellIS);
69: if (!cellIS) {
70: DMGetStratumIS(plex, "depth", depth, &cellIS);
71: }
72: DMGetLocalVector(plex, &locF);
73: VecZeroEntries(locF);
74: DMPlexComputeResidual_Internal(plex, cellIS, time, locX, NULL, time, locF, user);
75: DMLocalToGlobalBegin(plex, locF, ADD_VALUES, F);
76: DMLocalToGlobalEnd(plex, locF, ADD_VALUES, F);
77: DMRestoreLocalVector(plex, &locF);
78: ISDestroy(&cellIS);
79: DMDestroy(&plex);
80: return(0);
81: }
83: /*@
84: 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
86: Input Parameters:
87: + dm - The mesh
88: . t - The time
89: . locX - Local solution
90: . locX_t - Local solution time derivative, or NULL
91: - user - The user context
93: Level: developer
95: .seealso: DMPlexComputeJacobianActionFEM()
96: @*/
97: PetscErrorCode DMPlexTSComputeBoundary(DM dm, PetscReal time, Vec locX, Vec locX_t, void *user)
98: {
99: DM plex;
100: Vec faceGeometryFVM = NULL;
101: PetscInt Nf, f;
105: DMTSConvertPlex(dm, &plex, PETSC_TRUE);
106: DMGetNumFields(plex, &Nf);
107: if (!locX_t) {
108: /* This is the RHS part */
109: for (f = 0; f < Nf; f++) {
110: PetscObject obj;
111: PetscClassId id;
113: DMGetField(plex, f, NULL, &obj);
114: PetscObjectGetClassId(obj, &id);
115: if (id == PETSCFV_CLASSID) {
116: DMPlexGetGeometryFVM(plex, &faceGeometryFVM, NULL, NULL);
117: break;
118: }
119: }
120: }
121: DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, time, faceGeometryFVM, NULL, NULL);
122: DMPlexInsertTimeDerivativeBoundaryValues(plex, PETSC_TRUE, locX_t, time, faceGeometryFVM, NULL, NULL);
123: DMDestroy(&plex);
124: return(0);
125: }
127: /*@
128: DMPlexTSComputeIFunctionFEM - Form the local residual F from the local input X using pointwise functions specified by the user
130: Input Parameters:
131: + dm - The mesh
132: . t - The time
133: . locX - Local solution
134: . locX_t - Local solution time derivative, or NULL
135: - user - The user context
137: Output Parameter:
138: . locF - Local output vector
140: Level: developer
142: .seealso: DMPlexComputeJacobianActionFEM()
143: @*/
144: PetscErrorCode DMPlexTSComputeIFunctionFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, Vec locF, void *user)
145: {
146: DM plex;
147: IS cellIS;
148: PetscInt depth;
152: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
153: DMPlexGetDepth(plex, &depth);
154: DMGetStratumIS(plex, "dim", depth, &cellIS);
155: if (!cellIS) {
156: DMGetStratumIS(plex, "depth", depth, &cellIS);
157: }
158: DMPlexComputeResidual_Internal(plex, cellIS, time, locX, locX_t, time, locF, user);
159: ISDestroy(&cellIS);
160: DMDestroy(&plex);
161: return(0);
162: }
164: /*@
165: DMPlexTSComputeIJacobianFEM - Form the local Jacobian J from the local input X using pointwise functions specified by the user
167: Input Parameters:
168: + dm - The mesh
169: . t - The time
170: . locX - Local solution
171: . locX_t - Local solution time derivative, or NULL
172: . X_tshift - The multiplicative parameter for dF/du_t
173: - user - The user context
175: Output Parameter:
176: . locF - Local output vector
178: Level: developer
180: .seealso: DMPlexComputeJacobianActionFEM()
181: @*/
182: PetscErrorCode DMPlexTSComputeIJacobianFEM(DM dm, PetscReal time, Vec locX, Vec locX_t, PetscReal X_tShift, Mat Jac, Mat JacP, void *user)
183: {
184: DM plex;
185: PetscDS prob;
186: PetscBool hasJac, hasPrec;
187: IS cellIS;
188: PetscInt depth;
192: DMTSConvertPlex(dm,&plex,PETSC_TRUE);
193: DMGetDS(dm, &prob);
194: PetscDSHasJacobian(prob, &hasJac);
195: PetscDSHasJacobianPreconditioner(prob, &hasPrec);
196: if (hasJac && hasPrec) {MatZeroEntries(Jac);}
197: MatZeroEntries(JacP);
198: DMPlexGetDepth(plex,&depth);
199: DMGetStratumIS(plex, "dim", depth, &cellIS);
200: if (!cellIS) {DMGetStratumIS(plex, "depth", depth, &cellIS);}
201: DMPlexComputeJacobian_Internal(plex, cellIS, time, X_tShift, locX, locX_t, Jac, JacP, user);
202: ISDestroy(&cellIS);
203: DMDestroy(&plex);
204: return(0);
205: }
207: /*@C
208: DMTSCheckResidual - Check the residual of the exact solution
210: Input Parameters:
211: + ts - the TS object
212: . dm - the DM
213: . t - the time
214: . u - a DM vector
215: . u_t - a DM vector
216: - tol - A tolerance for the check, or -1 to print the results instead
218: Output Parameters:
219: . residual - The residual norm of the exact solution, or NULL
221: Level: developer
223: .seealso: DNTSCheckFromOptions(), DMTSCheckJacobian(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckJacobian()
224: @*/
225: PetscErrorCode DMTSCheckResidual(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscReal *residual)
226: {
227: MPI_Comm comm;
228: Vec r;
229: PetscReal res;
237: PetscObjectGetComm((PetscObject) ts, &comm);
238: DMComputeExactSolution(dm, t, u, u_t);
239: VecDuplicate(u, &r);
240: TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
241: VecNorm(r, NORM_2, &res);
242: if (tol >= 0.0) {
243: if (res > tol) SETERRQ2(comm, PETSC_ERR_ARG_WRONG, "L_2 Residual %g exceeds tolerance %g", (double) res, (double) tol);
244: } else if (residual) {
245: *residual = res;
246: } else {
247: PetscPrintf(comm, "L_2 Residual: %g\n", (double)res);
248: VecChop(r, 1.0e-10);
249: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", (PetscObject) dm);
250: PetscObjectSetName((PetscObject) r, "Initial Residual");
251: PetscObjectSetOptionsPrefix((PetscObject)r,"res_");
252: VecViewFromOptions(r, NULL, "-vec_view");
253: PetscObjectCompose((PetscObject) r, "__Vec_bc_zero__", NULL);
254: }
255: VecDestroy(&r);
256: return(0);
257: }
259: /*@C
260: DMTSCheckJacobian - Check the Jacobian of the exact solution against the residual using the Taylor Test
262: Input Parameters:
263: + ts - the TS object
264: . dm - the DM
265: . t - the time
266: . u - a DM vector
267: . u_t - a DM vector
268: - tol - A tolerance for the check, or -1 to print the results instead
270: Output Parameters:
271: + isLinear - Flag indicaing that the function looks linear, or NULL
272: - convRate - The rate of convergence of the linear model, or NULL
274: Level: developer
276: .seealso: DNTSCheckFromOptions(), DMTSCheckResidual(), DNSNESCheckFromOptions(), DMSNESCheckDiscretization(), DMSNESCheckResidual()
277: @*/
278: PetscErrorCode DMTSCheckJacobian(TS ts, DM dm, PetscReal t, Vec u, Vec u_t, PetscReal tol, PetscBool *isLinear, PetscReal *convRate)
279: {
280: MPI_Comm comm;
281: PetscDS ds;
282: Mat J, M;
283: MatNullSpace nullspace;
284: PetscReal dt, shift, slope, intercept;
285: PetscBool hasJac, hasPrec, isLin = PETSC_FALSE;
294: PetscObjectGetComm((PetscObject) ts, &comm);
295: DMComputeExactSolution(dm, t, u, u_t);
296: /* Create and view matrices */
297: TSGetTimeStep(ts, &dt);
298: shift = 1.0/dt;
299: DMCreateMatrix(dm, &J);
300: DMGetDS(dm, &ds);
301: PetscDSHasJacobian(ds, &hasJac);
302: PetscDSHasJacobianPreconditioner(ds, &hasPrec);
303: if (hasJac && hasPrec) {
304: DMCreateMatrix(dm, &M);
305: TSComputeIJacobian(ts, t, u, u_t, shift, J, M, PETSC_FALSE);
306: PetscObjectSetName((PetscObject) M, "Preconditioning Matrix");
307: PetscObjectSetOptionsPrefix((PetscObject) M, "jacpre_");
308: MatViewFromOptions(M, NULL, "-mat_view");
309: MatDestroy(&M);
310: } else {
311: TSComputeIJacobian(ts, t, u, u_t, shift, J, J, PETSC_FALSE);
312: }
313: PetscObjectSetName((PetscObject) J, "Jacobian");
314: PetscObjectSetOptionsPrefix((PetscObject) J, "jac_");
315: MatViewFromOptions(J, NULL, "-mat_view");
316: /* Check nullspace */
317: MatGetNullSpace(J, &nullspace);
318: if (nullspace) {
319: PetscBool isNull;
320: MatNullSpaceTest(nullspace, J, &isNull);
321: if (!isNull) SETERRQ(comm, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
322: }
323: MatNullSpaceDestroy(&nullspace);
324: /* Taylor test */
325: {
326: PetscRandom rand;
327: Vec du, uhat, uhat_t, r, rhat, df;
328: PetscReal h;
329: PetscReal *es, *hs, *errors;
330: PetscReal hMax = 1.0, hMin = 1e-6, hMult = 0.1;
331: PetscInt Nv, v;
333: /* Choose a perturbation direction */
334: PetscRandomCreate(comm, &rand);
335: VecDuplicate(u, &du);
336: VecSetRandom(du, rand);
337: PetscRandomDestroy(&rand);
338: VecDuplicate(u, &df);
339: MatMult(J, du, df);
340: /* Evaluate residual at u, F(u), save in vector r */
341: VecDuplicate(u, &r);
342: TSComputeIFunction(ts, t, u, u_t, r, PETSC_FALSE);
343: /* Look at the convergence of our Taylor approximation as we approach u */
344: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv);
345: PetscCalloc3(Nv, &es, Nv, &hs, Nv, &errors);
346: VecDuplicate(u, &uhat);
347: VecDuplicate(u, &uhat_t);
348: VecDuplicate(u, &rhat);
349: for (h = hMax, Nv = 0; h >= hMin; h *= hMult, ++Nv) {
350: VecWAXPY(uhat, h, du, u);
351: VecWAXPY(uhat_t, h*shift, du, u_t);
352: /* 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 */
353: TSComputeIFunction(ts, t, uhat, uhat_t, rhat, PETSC_FALSE);
354: VecAXPBYPCZ(rhat, -1.0, -h, 1.0, r, df);
355: VecNorm(rhat, NORM_2, &errors[Nv]);
357: es[Nv] = PetscLog10Real(errors[Nv]);
358: hs[Nv] = PetscLog10Real(h);
359: }
360: VecDestroy(&uhat);
361: VecDestroy(&uhat_t);
362: VecDestroy(&rhat);
363: VecDestroy(&df);
364: VecDestroy(&r);
365: VecDestroy(&du);
366: for (v = 0; v < Nv; ++v) {
367: if ((tol >= 0) && (errors[v] > tol)) break;
368: else if (errors[v] > PETSC_SMALL) break;
369: }
370: if (v == Nv) isLin = PETSC_TRUE;
371: PetscLinearRegression(Nv, hs, es, &slope, &intercept);
372: PetscFree3(es, hs, errors);
373: /* Slope should be about 2 */
374: if (tol >= 0) {
375: if (!isLin && PetscAbsReal(2 - slope) > tol) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Taylor approximation convergence rate should be 2, not %0.2f", (double) slope);
376: } else if (isLinear || convRate) {
377: if (isLinear) *isLinear = isLin;
378: if (convRate) *convRate = slope;
379: } else {
380: if (!isLin) {PetscPrintf(comm, "Taylor approximation converging at order %3.2f\n", (double) slope);}
381: else {PetscPrintf(comm, "Function appears to be linear\n");}
382: }
383: }
384: MatDestroy(&J);
385: return(0);
386: }
388: /*@C
389: DMTSCheckFromOptions - Check the residual and Jacobian functions using the exact solution by outputting some diagnostic information
391: Input Parameters:
392: + ts - the TS object
393: - u - representative TS vector
395: Note: The user must call PetscDSSetExactSolution() beforehand
397: Level: developer
398: @*/
399: PetscErrorCode DMTSCheckFromOptions(TS ts, Vec u)
400: {
401: DM dm;
402: SNES snes;
403: Vec sol, u_t;
404: PetscReal t;
405: PetscBool check;
409: PetscOptionsHasName(((PetscObject)ts)->options,((PetscObject)ts)->prefix, "-dmts_check", &check);
410: if (!check) return(0);
411: VecDuplicate(u, &sol);
412: VecCopy(u, sol);
413: TSSetSolution(ts, u);
414: TSGetDM(ts, &dm);
415: TSSetUp(ts);
416: TSGetSNES(ts, &snes);
417: SNESSetSolution(snes, u);
419: TSGetTime(ts, &t);
420: DMSNESCheckDiscretization(snes, dm, t, sol, -1.0, NULL);
421: DMGetGlobalVector(dm, &u_t);
422: DMTSCheckResidual(ts, dm, t, sol, u_t, -1.0, NULL);
423: DMTSCheckJacobian(ts, dm, t, sol, u_t, -1.0, NULL, NULL);
424: DMRestoreGlobalVector(dm, &u_t);
426: VecDestroy(&sol);
427: return(0);
428: }