1 #include <petscconvest.h> /*I "petscconvest.h" I*/
2 #include <petscts.h>
3 #include <petscdmplex.h>
4
5 #include <petsc/private/petscconvestimpl.h>
6
PetscConvEstSetTS_Private(PetscConvEst ce,PetscObject solver)7 static PetscErrorCode PetscConvEstSetTS_Private(PetscConvEst ce, PetscObject solver)
8 {
9 PetscClassId id;
10
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 }
17
PetscConvEstInitGuessTS_Private(PetscConvEst ce,PetscInt r,DM dm,Vec u)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 }
24
PetscConvEstComputeErrorTS_Private(PetscConvEst ce,PetscInt r,DM dm,Vec u,PetscReal errors[])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);
29
30 PetscFunctionBegin;
31 PetscCall(TSGetComputeExactError(ts, &exactError));
32 if (exactError) {
33 Vec e;
34 PetscInt f;
35
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;
43
44 PetscCall(TSGetSolveTime(ts, &t));
45 PetscCall(DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors));
46 }
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce,PetscReal alpha[])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;
56
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 }
115
PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce,PetscReal alpha[])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;
128
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));
138
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;
146
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;
157
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 *);
169
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;
208
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;
224
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;
257
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 }
278
PetscConvEstUseTS(PetscConvEst ce,PetscBool checkTemporal)279 PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
280 {
281 PetscFunctionBegin;
282 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1);
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 }
293