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: void *fctx, *jctx, *rctx;
128: void *ctx;
130: PetscFunctionBegin;
131: PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r);
132: PetscCall(DMGetDimension(ce->idm, &dim));
133: PetscCall(DMGetApplicationContext(ce->idm, &ctx));
134: PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE));
135: PetscCall(DMGetRefineLevel(ce->idm, &oldlevel));
136: PetscCall(PetscMalloc1((Nr + 1), &dm));
137: PetscCall(TSGetSolution(ts, &uInitial));
138: /* Loop over meshes */
139: dm[0] = ce->idm;
140: for (r = 0; r <= Nr; ++r) {
141: Vec u;
142: #if defined(PETSC_USE_LOG)
143: PetscLogStage stage;
144: #endif
145: char stageName[PETSC_MAX_PATH_LEN];
146: const char *dmname, *uname;
148: PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r));
149: #if defined(PETSC_USE_LOG)
150: PetscCall(PetscLogStageGetId(stageName, &stage));
151: if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage));
152: #endif
153: PetscCall(PetscLogStagePush(stage));
154: if (r > 0) {
155: if (!ce->noRefine) {
156: PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]));
157: PetscCall(DMSetCoarseDM(dm[r], dm[r - 1]));
158: } else {
159: DM cdm, rcdm;
161: PetscCall(DMClone(dm[r - 1], &dm[r]));
162: PetscCall(DMCopyDisc(dm[r - 1], dm[r]));
163: PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm));
164: PetscCall(DMGetCoordinateDM(dm[r], &rcdm));
165: PetscCall(DMCopyDisc(cdm, rcdm));
166: }
167: PetscCall(DMCopyTransform(ce->idm, dm[r]));
168: PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname));
169: PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname));
170: for (f = 0; f <= Nf; ++f) {
171: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
173: PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr));
174: PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr));
175: }
176: }
177: PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view"));
178: /* Create solution */
179: PetscCall(DMCreateGlobalVector(dm[r], &u));
180: PetscCall(DMGetField(dm[r], 0, NULL, &disc));
181: PetscCall(PetscObjectGetName(disc, &uname));
182: PetscCall(PetscObjectSetName((PetscObject)u, uname));
183: /* Setup solver */
184: PetscCall(TSReset(ts));
185: PetscCall(TSSetDM(ts, dm[r]));
186: PetscCall(DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx));
187: if (r > 0) {
188: PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx));
189: PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx));
190: PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx));
191: if (ifunc) PetscCall(DMTSSetIFunctionLocal(dm[r], ifunc, fctx));
192: if (ijac) PetscCall(DMTSSetIJacobianLocal(dm[r], ijac, jctx));
193: if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(dm[r], rhsfunc, rctx));
194: }
195: PetscCall(TSSetTime(ts, 0.0));
196: PetscCall(TSSetStepNumber(ts, 0));
197: PetscCall(TSSetFromOptions(ts));
198: PetscCall(TSSetSolution(ts, u));
199: PetscCall(VecDestroy(&u));
200: /* Create initial guess */
201: PetscCall(TSGetSolution(ts, &u));
202: PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u));
203: PetscCall(TSSolve(ts, NULL));
204: PetscCall(TSGetSolution(ts, &u));
205: PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
206: PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * Nf]));
207: PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
208: for (f = 0; f < Nf; ++f) {
209: PetscSection s, fs;
210: PetscInt lsize;
212: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
213: PetscCall(DMGetLocalSection(dm[r], &s));
214: PetscCall(PetscSectionGetField(s, f, &fs));
215: PetscCall(PetscSectionGetConstrainedStorageSize(fs, &lsize));
216: PetscCall(MPIU_Allreduce(&lsize, &ce->dofs[r * Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)ts)));
217: PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]));
218: PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]));
219: }
220: /* Monitor */
221: PetscCall(PetscConvEstMonitorDefault(ce, r));
222: if (!r) {
223: /* PCReset() does not wipe out the level structure */
224: SNES snes;
225: KSP ksp;
226: PC pc;
228: PetscCall(TSGetSNES(ts, &snes));
229: PetscCall(SNESGetKSP(snes, &ksp));
230: PetscCall(KSPGetPC(ksp, &pc));
231: PetscCall(PCMGGetLevels(pc, &oldnlev));
232: }
233: /* Cleanup */
234: PetscCall(PetscLogStagePop());
235: }
236: PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx));
237: PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx));
238: PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx));
239: for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r]));
240: /* Fit convergence rate */
241: PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
242: for (f = 0; f < Nf; ++f) {
243: for (r = 0; r <= Nr; ++r) {
244: x[r] = PetscLog10Real(ce->dofs[r * Nf + f]);
245: y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
246: }
247: PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
248: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
249: alpha[f] = -slope * dim;
250: }
251: PetscCall(PetscFree2(x, y));
252: PetscCall(PetscFree(dm));
253: /* Restore solver */
254: PetscCall(TSReset(ts));
255: {
256: /* PCReset() does not wipe out the level structure */
257: SNES snes;
258: KSP ksp;
259: PC pc;
261: PetscCall(TSGetSNES(ts, &snes));
262: PetscCall(SNESGetKSP(snes, &ksp));
263: PetscCall(KSPGetPC(ksp, &pc));
264: PetscCall(PCMGSetLevels(pc, oldnlev, NULL));
265: PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */
266: }
267: PetscCall(TSSetDM(ts, ce->idm));
268: PetscCall(DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx));
269: if (ifunc) PetscCall(DMTSSetIFunctionLocal(ce->idm, ifunc, fctx));
270: if (ijac) PetscCall(DMTSSetIJacobianLocal(ce->idm, ijac, jctx));
271: if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(ce->idm, rhsfunc, rctx));
272: PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
273: PetscCall(TSSetTime(ts, 0.0));
274: PetscCall(TSSetStepNumber(ts, 0));
275: PetscCall(TSSetFromOptions(ts));
276: PetscCall(TSSetSolution(ts, uInitial));
277: PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial));
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
282: {
283: PetscFunctionBegin;
285: ce->ops->setsolver = PetscConvEstSetTS_Private;
286: ce->ops->initguess = PetscConvEstInitGuessTS_Private;
287: ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
288: if (checkTemporal) {
289: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
290: } else {
291: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
292: }
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }