Actual source code: tsconvest.c
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;
11: PetscFunctionBegin;
12: PetscCall(PetscObjectGetClassId(ce->solver, &id));
13: PetscCheck(id == TS_CLASSID, PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS");
14: PetscCall(TSGetDM((TS)ce->solver, &ce->idm));
15: PetscFunctionReturn(PETSC_SUCCESS);
16: }
18: static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
19: {
20: PetscFunctionBegin;
21: PetscCall(TSComputeInitialCondition((TS)ce->solver, u));
22: PetscFunctionReturn(PETSC_SUCCESS);
23: }
25: static PetscErrorCode PetscConvEstComputeErrorTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
26: {
27: TS ts = (TS)ce->solver;
28: PetscErrorCode (*exactError)(TS, Vec, Vec);
30: PetscFunctionBegin;
31: PetscCall(TSGetComputeExactError(ts, &exactError));
32: if (exactError) {
33: Vec e;
34: PetscInt f;
36: PetscCall(VecDuplicate(u, &e));
37: PetscCall(TSComputeExactError(ts, u, e));
38: PetscCall(VecNorm(e, NORM_2, errors));
39: for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
40: PetscCall(VecDestroy(&e));
41: } else {
42: PetscReal t;
44: PetscCall(TSGetSolveTime(ts, &t));
45: PetscCall(DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors));
46: }
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce, PetscReal alpha[])
51: {
52: TS ts = (TS)ce->solver;
53: Vec u, u0;
54: PetscReal *dt, *x, *y, slope, intercept;
55: PetscInt Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r;
57: PetscFunctionBegin;
58: PetscCall(PetscMalloc1(Nr + 1, &dt));
59: PetscCall(TSGetTimeStep(ts, &dt[0]));
60: PetscCall(TSGetMaxSteps(ts, &oNs));
61: PetscCall(TSGetSolution(ts, &u0));
62: PetscCall(PetscObjectReference((PetscObject)u0));
63: Ns = oNs;
64: for (r = 0; r <= Nr; ++r) {
65: if (r > 0) {
66: dt[r] = dt[r - 1] / ce->r;
67: Ns = PetscCeilReal(Ns * ce->r);
68: }
69: PetscCall(TSSetTime(ts, 0.0));
70: PetscCall(TSSetStepNumber(ts, 0));
71: PetscCall(TSSetTimeStep(ts, dt[r]));
72: PetscCall(TSSetMaxSteps(ts, Ns));
73: PetscCall(TSGetSolution(ts, &u));
74: PetscCall(PetscConvEstComputeInitialGuess(ce, r, NULL, u));
75: PetscCall(TSSolve(ts, NULL));
76: PetscCall(TSGetSolution(ts, &u));
77: PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
78: PetscCall(PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r * Nf]));
79: PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
80: for (f = 0; f < Nf; ++f) {
81: ce->dofs[r * Nf + f] = 1.0 / dt[r];
82: PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]));
83: PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]));
84: }
85: /* Monitor */
86: PetscCall(PetscConvEstMonitorDefault(ce, r));
87: }
88: /* Fit convergence rate */
89: if (Nr) {
90: PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
91: for (f = 0; f < Nf; ++f) {
92: for (r = 0; r <= Nr; ++r) {
93: x[r] = PetscLog10Real(dt[r]);
94: y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
95: }
96: PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
97: /* Since lg err = s lg dt + b */
98: alpha[f] = slope;
99: }
100: PetscCall(PetscFree2(x, y));
101: }
102: /* Reset solver */
103: PetscCall(TSReset(ts));
104: PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
105: PetscCall(TSSetTime(ts, 0.0));
106: PetscCall(TSSetStepNumber(ts, 0));
107: PetscCall(TSSetTimeStep(ts, dt[0]));
108: PetscCall(TSSetMaxSteps(ts, oNs));
109: PetscCall(TSSetSolution(ts, u0));
110: PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, u0));
111: PetscCall(VecDestroy(&u0));
112: PetscCall(PetscFree(dt));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce, PetscReal alpha[])
117: {
118: TS ts = (TS)ce->solver;
119: Vec uInitial;
120: DM *dm;
121: PetscObject disc;
122: PetscReal *x, *y, slope, intercept;
123: PetscInt Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev;
124: PetscErrorCode (*ifunc)(DM, PetscReal, Vec, Vec, Vec, void *);
125: PetscErrorCode (*ijac)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
126: PetscErrorCode (*rhsfunc)(DM, PetscReal, Vec, Vec, void *);
127: PetscCtx fctx, jctx, rctx, ctx;
129: PetscFunctionBegin;
130: PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r);
131: PetscCall(DMGetDimension(ce->idm, &dim));
132: PetscCall(DMGetApplicationContext(ce->idm, &ctx));
133: PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE));
134: PetscCall(DMGetRefineLevel(ce->idm, &oldlevel));
135: PetscCall(PetscMalloc1(Nr + 1, &dm));
136: PetscCall(TSGetSolution(ts, &uInitial));
137: PetscCall(PetscObjectReference((PetscObject)uInitial));
139: /* Loop over meshes */
140: dm[0] = ce->idm;
141: for (r = 0; r <= Nr; ++r) {
142: Vec u;
143: PetscLogStage stage;
144: char stageName[PETSC_MAX_PATH_LEN];
145: const char *dmname, *uname;
147: PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r));
148: PetscCall(PetscLogStageGetId(stageName, &stage));
149: if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage));
150: PetscCall(PetscLogStagePush(stage));
151: if (r > 0) {
152: if (!ce->noRefine) {
153: PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]));
154: PetscCall(DMSetCoarseDM(dm[r], dm[r - 1]));
155: } else {
156: DM cdm, rcdm;
158: PetscCall(DMClone(dm[r - 1], &dm[r]));
159: PetscCall(DMCopyDisc(dm[r - 1], dm[r]));
160: PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm));
161: PetscCall(DMGetCoordinateDM(dm[r], &rcdm));
162: PetscCall(DMCopyDisc(cdm, rcdm));
163: }
164: PetscCall(DMCopyTransform(ce->idm, dm[r]));
165: PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname));
166: PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname));
167: for (f = 0; f < Nf; ++f) {
168: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
170: PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr));
171: PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr));
172: }
173: }
174: PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view"));
175: /* Create solution */
176: PetscCall(DMCreateGlobalVector(dm[r], &u));
177: PetscCall(DMGetField(dm[r], 0, NULL, &disc));
178: PetscCall(PetscObjectGetName(disc, &uname));
179: PetscCall(PetscObjectSetName((PetscObject)u, uname));
180: /* Setup solver */
181: PetscCall(TSReset(ts));
182: PetscCall(TSSetDM(ts, dm[r]));
183: PetscCall(DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx));
184: if (r > 0) {
185: PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx));
186: PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx));
187: PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx));
188: if (ifunc) PetscCall(DMTSSetIFunctionLocal(dm[r], ifunc, fctx));
189: if (ijac) PetscCall(DMTSSetIJacobianLocal(dm[r], ijac, jctx));
190: if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(dm[r], rhsfunc, rctx));
191: }
192: PetscCall(TSSetTime(ts, 0.0));
193: PetscCall(TSSetStepNumber(ts, 0));
194: PetscCall(TSSetFromOptions(ts));
195: PetscCall(TSSetSolution(ts, u));
196: PetscCall(VecDestroy(&u));
197: /* Create initial guess */
198: PetscCall(TSGetSolution(ts, &u));
199: PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u));
200: PetscCall(TSSolve(ts, NULL));
201: PetscCall(TSGetSolution(ts, &u));
202: PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
203: PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * Nf]));
204: PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
205: for (f = 0; f < Nf; ++f) {
206: PetscSection s, fs;
207: PetscInt lsize;
209: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
210: PetscCall(DMGetLocalSection(dm[r], &s));
211: PetscCall(PetscSectionGetField(s, f, &fs));
212: PetscCall(PetscSectionGetConstrainedStorageSize(fs, &lsize));
213: PetscCallMPI(MPIU_Allreduce(&lsize, &ce->dofs[r * Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)ts)));
214: PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]));
215: PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]));
216: }
217: /* Monitor */
218: PetscCall(PetscConvEstMonitorDefault(ce, r));
219: if (!r) {
220: /* PCReset() does not wipe out the level structure */
221: SNES snes;
222: KSP ksp;
223: PC pc;
225: PetscCall(TSGetSNES(ts, &snes));
226: PetscCall(SNESGetKSP(snes, &ksp));
227: PetscCall(KSPGetPC(ksp, &pc));
228: PetscCall(PCMGGetLevels(pc, &oldnlev));
229: }
230: /* Cleanup */
231: PetscCall(PetscLogStagePop());
232: }
233: PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx));
234: PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx));
235: PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx));
236: for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r]));
237: /* Fit convergence rate */
238: PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
239: for (f = 0; f < Nf; ++f) {
240: for (r = 0; r <= Nr; ++r) {
241: x[r] = PetscLog10Real(ce->dofs[r * Nf + f]);
242: y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
243: }
244: PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
245: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
246: alpha[f] = -slope * dim;
247: }
248: PetscCall(PetscFree2(x, y));
249: PetscCall(PetscFree(dm));
250: /* Restore solver */
251: PetscCall(TSReset(ts));
252: {
253: /* PCReset() does not wipe out the level structure */
254: SNES snes;
255: KSP ksp;
256: PC pc;
258: PetscCall(TSGetSNES(ts, &snes));
259: PetscCall(SNESGetKSP(snes, &ksp));
260: PetscCall(KSPGetPC(ksp, &pc));
261: PetscCall(PCMGSetLevels(pc, oldnlev, NULL));
262: PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */
263: }
264: PetscCall(TSSetDM(ts, ce->idm));
265: PetscCall(DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx));
266: if (ifunc) PetscCall(DMTSSetIFunctionLocal(ce->idm, ifunc, fctx));
267: if (ijac) PetscCall(DMTSSetIJacobianLocal(ce->idm, ijac, jctx));
268: if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(ce->idm, rhsfunc, rctx));
269: PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
270: PetscCall(TSSetTime(ts, 0.0));
271: PetscCall(TSSetStepNumber(ts, 0));
272: PetscCall(TSSetFromOptions(ts));
273: PetscCall(TSSetSolution(ts, uInitial));
274: PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial));
275: PetscCall(VecDestroy(&uInitial));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
280: {
281: PetscFunctionBegin;
283: ce->ops->setsolver = PetscConvEstSetTS_Private;
284: ce->ops->initguess = PetscConvEstInitGuessTS_Private;
285: ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
286: if (checkTemporal) {
287: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
288: } else {
289: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
290: }
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }