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