Actual source code: tsconvest.c
petsc-3.14.6 2021-03-30
1: #include <petscconvest.h>
2: #include <petscts.h>
3: #include <petscdmplex.h>
5: #include <petsc/private/petscconvestimpl.h>
7: static PetscErrorCode PetscConvEstSetTS_Private(PetscConvEst ce, PetscObject solver)
8: {
9: PetscClassId id;
13: PetscObjectGetClassId(ce->solver, &id);
14: if (id != TS_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS");
15: TSGetDM((TS) ce->solver, &ce->idm);
16: return(0);
17: }
19: static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
20: {
24: TSComputeInitialCondition((TS) ce->solver, u);
25: return(0);
26: }
28: static PetscErrorCode PetscConvEstComputeErrorTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
29: {
30: TS ts = (TS) ce->solver;
31: PetscErrorCode (*exactError)(TS, Vec, Vec);
32: PetscErrorCode ierr;
35: TSGetComputeExactError(ts, &exactError);
36: if (exactError) {
37: Vec e;
38: PetscInt f;
40: VecDuplicate(u, &e);
41: TSComputeExactError(ts, u, e);
42: VecNorm(e, NORM_2, errors);
43: for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
44: VecDestroy(&e);
45: } else {
46: PetscReal t;
48: TSGetSolveTime(ts, &t);
49: DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors);
50: }
51: return(0);
52: }
54: static PetscErrorCode PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce, PetscReal alpha[])
55: {
56: TS ts = (TS) ce->solver;
57: Vec u;
58: PetscReal *dt, *x, *y, slope, intercept;
59: PetscInt Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r;
63: TSGetSolution(ts, &u);
64: PetscMalloc1(Nr+1, &dt);
65: TSGetTimeStep(ts, &dt[0]);
66: TSGetMaxSteps(ts, &oNs);
67: Ns = oNs;
68: for (r = 0; r <= Nr; ++r) {
69: if (r > 0) {
70: dt[r] = dt[r-1]/ce->r;
71: Ns = PetscCeilReal(Ns*ce->r);
72: }
73: TSSetTime(ts, 0.0);
74: TSSetStepNumber(ts, 0);
75: TSSetTimeStep(ts, dt[r]);
76: TSSetMaxSteps(ts, Ns);
77: PetscConvEstComputeInitialGuess(ce, r, NULL, u);
78: TSSolve(ts, u);
79: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
80: PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r*Nf]);
81: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
82: for (f = 0; f < Nf; ++f) {
83: PetscLogEventSetDof(ce->event, f, 1.0/dt[r]);
84: PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
85: }
86: /* Monitor */
87: PetscConvEstMonitorDefault(ce, r);
88: }
89: /* Fit convergence rate */
90: if (Nr) {
91: PetscMalloc2(Nr+1, &x, Nr+1, &y);
92: for (f = 0; f < Nf; ++f) {
93: for (r = 0; r <= Nr; ++r) {
94: x[r] = PetscLog10Real(dt[r]);
95: y[r] = PetscLog10Real(ce->errors[r*Nf+f]);
96: }
97: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
98: /* Since lg err = s lg dt + b */
99: alpha[f] = slope;
100: }
101: PetscFree2(x, y);
102: }
103: /* Reset solver */
104: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
105: TSSetTime(ts, 0.0);
106: TSSetStepNumber(ts, 0);
107: TSSetTimeStep(ts, dt[0]);
108: TSSetMaxSteps(ts, oNs);
109: PetscConvEstComputeInitialGuess(ce, 0, NULL, u);
110: PetscFree(dt);
111: return(0);
112: }
114: static PetscErrorCode PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce, PetscReal alpha[])
115: {
116: TS ts = (TS) ce->solver;
117: Vec uInitial;
118: DM *dm;
119: PetscObject disc;
120: PetscReal *x, *y, slope, intercept;
121: PetscInt *dof, Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev;
122: void *ctx;
126: if (ce->r != 2.0) SETERRQ1(PetscObjectComm((PetscObject) ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double) ce->r);
127: DMGetDimension(ce->idm, &dim);
128: DMGetApplicationContext(ce->idm, &ctx);
129: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
130: DMGetRefineLevel(ce->idm, &oldlevel);
131: PetscMalloc2((Nr+1), &dm, (Nr+1)*Nf, &dof);
132: TSGetSolution(ts, &uInitial);
133: /* Loop over meshes */
134: dm[0] = ce->idm;
135: for (r = 0; r <= Nr; ++r) {
136: Vec u;
137: #if defined(PETSC_USE_LOG)
138: PetscLogStage stage;
139: #endif
140: char stageName[PETSC_MAX_PATH_LEN];
141: const char *dmname, *uname;
143: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
144: PetscLogStageRegister(stageName, &stage);
145: PetscLogStagePush(stage);
146: if (r > 0) {
147: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
148: DMSetCoarseDM(dm[r], dm[r-1]);
149: DMCopyTransform(ce->idm, dm[r]);
150: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
151: PetscObjectSetName((PetscObject) dm[r], dmname);
152: for (f = 0; f <= Nf; ++f) {
153: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
155: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
156: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
157: }
158: }
159: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
160: /* Create solution */
161: DMCreateGlobalVector(dm[r], &u);
162: DMGetField(dm[r], 0, NULL, &disc);
163: PetscObjectGetName(disc, &uname);
164: PetscObjectSetName((PetscObject) u, uname);
165: /* Setup solver */
166: TSReset(ts);
167: TSSetDM(ts, dm[r]);
168: DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx);
169: DMTSSetIFunctionLocal(dm[r], DMPlexTSComputeIFunctionFEM, ctx);
170: DMTSSetIJacobianLocal(dm[r], DMPlexTSComputeIJacobianFEM, ctx);
171: TSSetTime(ts, 0.0);
172: TSSetStepNumber(ts, 0);
173: TSSetFromOptions(ts);
174: /* Create initial guess */
175: PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
176: TSSolve(ts, u);
177: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
178: PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*Nf]);
179: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
180: for (f = 0; f < Nf; ++f) {
181: PetscSection s, fs;
182: PetscInt lsize;
184: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
185: DMGetLocalSection(dm[r], &s);
186: PetscSectionGetField(s, f, &fs);
187: PetscSectionGetConstrainedStorageSize(fs, &lsize);
188: MPI_Allreduce(&lsize, &dof[r*Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ts));
189: PetscLogEventSetDof(ce->event, f, dof[r*Nf+f]);
190: PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
191: }
192: /* Monitor */
193: PetscConvEstMonitorDefault(ce, r);
194: if (!r) {
195: /* PCReset() does not wipe out the level structure */
196: SNES snes;
197: KSP ksp;
198: PC pc;
200: TSGetSNES(ts, &snes);
201: SNESGetKSP(snes, &ksp);
202: KSPGetPC(ksp, &pc);
203: PCMGGetLevels(pc, &oldnlev);
204: }
205: /* Cleanup */
206: VecDestroy(&u);
207: PetscLogStagePop();
208: }
209: for (r = 1; r <= Nr; ++r) {
210: DMDestroy(&dm[r]);
211: }
212: /* Fit convergence rate */
213: PetscMalloc2(Nr+1, &x, Nr+1, &y);
214: for (f = 0; f < Nf; ++f) {
215: for (r = 0; r <= Nr; ++r) {
216: x[r] = PetscLog10Real(dof[r*Nf+f]);
217: y[r] = PetscLog10Real(ce->errors[r*Nf+f]);
218: }
219: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
220: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
221: alpha[f] = -slope * dim;
222: }
223: PetscFree2(x, y);
224: PetscFree2(dm, dof);
225: /* Restore solver */
226: TSReset(ts);
227: {
228: /* PCReset() does not wipe out the level structure */
229: SNES snes;
230: KSP ksp;
231: PC pc;
233: TSGetSNES(ts, &snes);
234: SNESGetKSP(snes, &ksp);
235: KSPGetPC(ksp, &pc);
236: PCMGSetLevels(pc, oldnlev, NULL);
237: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
238: }
239: TSSetDM(ts, ce->idm);
240: DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx);
241: DMTSSetIFunctionLocal(ce->idm, DMPlexTSComputeIFunctionFEM, ctx);
242: DMTSSetIJacobianLocal(ce->idm, DMPlexTSComputeIJacobianFEM, ctx);
243: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
244: TSSetTime(ts, 0.0);
245: TSSetStepNumber(ts, 0);
246: TSSetFromOptions(ts);
247: TSSetSolution(ts, uInitial);
248: PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial);
249: return(0);
250: }
252: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
253: {
256: ce->ops->setsolver = PetscConvEstSetTS_Private;
257: ce->ops->initguess = PetscConvEstInitGuessTS_Private;
258: ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
259: if (checkTemporal) {
260: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
261: } else {
262: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
263: }
264: return(0);
265: }