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 PetscErrorCode ierr;
11
12 PetscFunctionBegin;
13 ierr = PetscObjectGetClassId(ce->solver, &id);CHKERRQ(ierr);
14 if (id != TS_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS");
15 ierr = TSGetDM((TS) ce->solver, &ce->idm);CHKERRQ(ierr);
16 PetscFunctionReturn(0);
17 }
18
PetscConvEstInitGuessTS_Private(PetscConvEst ce,PetscInt r,DM dm,Vec u)19 static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
20 {
21 PetscErrorCode ierr;
22
23 PetscFunctionBegin;
24 ierr = TSComputeInitialCondition((TS) ce->solver, u);CHKERRQ(ierr);
25 PetscFunctionReturn(0);
26 }
27
PetscConvEstComputeErrorTS_Private(PetscConvEst ce,PetscInt r,DM dm,Vec u,PetscReal errors[])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;
33
34 PetscFunctionBegin;
35 ierr = TSGetComputeExactError(ts, &exactError);CHKERRQ(ierr);
36 if (exactError) {
37 Vec e;
38 PetscInt f;
39
40 ierr = VecDuplicate(u, &e);CHKERRQ(ierr);
41 ierr = TSComputeExactError(ts, u, e);CHKERRQ(ierr);
42 ierr = VecNorm(e, NORM_2, errors);CHKERRQ(ierr);
43 for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
44 ierr = VecDestroy(&e);CHKERRQ(ierr);
45 } else {
46 PetscReal t;
47
48 ierr = TSGetSolveTime(ts, &t);CHKERRQ(ierr);
49 ierr = DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors);CHKERRQ(ierr);
50 }
51 PetscFunctionReturn(0);
52 }
53
PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce,PetscReal alpha[])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;
60 PetscErrorCode ierr;
61
62 PetscFunctionBegin;
63 ierr = TSGetSolution(ts, &u);CHKERRQ(ierr);
64 ierr = PetscMalloc1(Nr+1, &dt);CHKERRQ(ierr);
65 ierr = TSGetTimeStep(ts, &dt[0]);CHKERRQ(ierr);
66 ierr = TSGetMaxSteps(ts, &oNs);CHKERRQ(ierr);
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 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr);
74 ierr = TSSetStepNumber(ts, 0);CHKERRQ(ierr);
75 ierr = TSSetTimeStep(ts, dt[r]);CHKERRQ(ierr);
76 ierr = TSSetMaxSteps(ts, Ns);CHKERRQ(ierr);
77 ierr = PetscConvEstComputeInitialGuess(ce, r, NULL, u);CHKERRQ(ierr);
78 ierr = TSSolve(ts, u);CHKERRQ(ierr);
79 ierr = PetscLogEventBegin(ce->event, ce, 0, 0, 0);CHKERRQ(ierr);
80 ierr = PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r*Nf]);CHKERRQ(ierr);
81 ierr = PetscLogEventEnd(ce->event, ce, 0, 0, 0);CHKERRQ(ierr);
82 for (f = 0; f < Nf; ++f) {
83 ierr = PetscLogEventSetDof(ce->event, f, 1.0/dt[r]);CHKERRQ(ierr);
84 ierr = PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);CHKERRQ(ierr);
85 }
86 /* Monitor */
87 ierr = PetscConvEstMonitorDefault(ce, r);CHKERRQ(ierr);
88 }
89 /* Fit convergence rate */
90 if (Nr) {
91 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr);
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 ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr);
98 /* Since lg err = s lg dt + b */
99 alpha[f] = slope;
100 }
101 ierr = PetscFree2(x, y);CHKERRQ(ierr);
102 }
103 /* Reset solver */
104 ierr = TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);CHKERRQ(ierr);
105 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr);
106 ierr = TSSetStepNumber(ts, 0);CHKERRQ(ierr);
107 ierr = TSSetTimeStep(ts, dt[0]);CHKERRQ(ierr);
108 ierr = TSSetMaxSteps(ts, oNs);CHKERRQ(ierr);
109 ierr = PetscConvEstComputeInitialGuess(ce, 0, NULL, u);CHKERRQ(ierr);
110 ierr = PetscFree(dt);CHKERRQ(ierr);
111 PetscFunctionReturn(0);
112 }
113
PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce,PetscReal alpha[])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;
123 PetscErrorCode ierr;
124
125 PetscFunctionBegin;
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 ierr = DMGetDimension(ce->idm, &dim);CHKERRQ(ierr);
128 ierr = DMGetApplicationContext(ce->idm, &ctx);CHKERRQ(ierr);
129 ierr = DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);CHKERRQ(ierr);
130 ierr = DMGetRefineLevel(ce->idm, &oldlevel);CHKERRQ(ierr);
131 ierr = PetscMalloc2((Nr+1), &dm, (Nr+1)*Nf, &dof);CHKERRQ(ierr);
132 ierr = TSGetSolution(ts, &uInitial);CHKERRQ(ierr);
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;
142
143 ierr = PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);CHKERRQ(ierr);
144 ierr = PetscLogStageRegister(stageName, &stage);CHKERRQ(ierr);
145 ierr = PetscLogStagePush(stage);CHKERRQ(ierr);
146 if (r > 0) {
147 ierr = DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);CHKERRQ(ierr);
148 ierr = DMSetCoarseDM(dm[r], dm[r-1]);CHKERRQ(ierr);
149 ierr = DMCopyTransform(ce->idm, dm[r]);CHKERRQ(ierr);
150 ierr = PetscObjectGetName((PetscObject) dm[r-1], &dmname);CHKERRQ(ierr);
151 ierr = PetscObjectSetName((PetscObject) dm[r], dmname);CHKERRQ(ierr);
152 for (f = 0; f <= Nf; ++f) {
153 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
154
155 ierr = DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);CHKERRQ(ierr);
156 ierr = DMSetNullSpaceConstructor(dm[r], f, nspconstr);CHKERRQ(ierr);
157 }
158 }
159 ierr = DMViewFromOptions(dm[r], NULL, "-conv_dm_view");CHKERRQ(ierr);
160 /* Create solution */
161 ierr = DMCreateGlobalVector(dm[r], &u);CHKERRQ(ierr);
162 ierr = DMGetField(dm[r], 0, NULL, &disc);CHKERRQ(ierr);
163 ierr = PetscObjectGetName(disc, &uname);CHKERRQ(ierr);
164 ierr = PetscObjectSetName((PetscObject) u, uname);CHKERRQ(ierr);
165 /* Setup solver */
166 ierr = TSReset(ts);CHKERRQ(ierr);
167 ierr = TSSetDM(ts, dm[r]);CHKERRQ(ierr);
168 ierr = DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx);CHKERRQ(ierr);
169 ierr = DMTSSetIFunctionLocal(dm[r], DMPlexTSComputeIFunctionFEM, ctx);CHKERRQ(ierr);
170 ierr = DMTSSetIJacobianLocal(dm[r], DMPlexTSComputeIJacobianFEM, ctx);CHKERRQ(ierr);
171 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr);
172 ierr = TSSetStepNumber(ts, 0);CHKERRQ(ierr);
173 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
174 /* Create initial guess */
175 ierr = PetscConvEstComputeInitialGuess(ce, r, dm[r], u);CHKERRQ(ierr);
176 ierr = TSSolve(ts, u);CHKERRQ(ierr);
177 ierr = PetscLogEventBegin(ce->event, ce, 0, 0, 0);CHKERRQ(ierr);
178 ierr = PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*Nf]);CHKERRQ(ierr);
179 ierr = PetscLogEventEnd(ce->event, ce, 0, 0, 0);CHKERRQ(ierr);
180 for (f = 0; f < Nf; ++f) {
181 PetscSection s, fs;
182 PetscInt lsize;
183
184 /* Could use DMGetOutputDM() to add in Dirichlet dofs */
185 ierr = DMGetLocalSection(dm[r], &s);CHKERRQ(ierr);
186 ierr = PetscSectionGetField(s, f, &fs);CHKERRQ(ierr);
187 ierr = PetscSectionGetConstrainedStorageSize(fs, &lsize);CHKERRQ(ierr);
188 ierr = MPI_Allreduce(&lsize, &dof[r*Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ts));CHKERRQ(ierr);
189 ierr = PetscLogEventSetDof(ce->event, f, dof[r*Nf+f]);CHKERRQ(ierr);
190 ierr = PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);CHKERRQ(ierr);
191 }
192 /* Monitor */
193 ierr = PetscConvEstMonitorDefault(ce, r);CHKERRQ(ierr);
194 if (!r) {
195 /* PCReset() does not wipe out the level structure */
196 SNES snes;
197 KSP ksp;
198 PC pc;
199
200 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr);
201 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
202 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
203 ierr = PCMGGetLevels(pc, &oldnlev);CHKERRQ(ierr);
204 }
205 /* Cleanup */
206 ierr = VecDestroy(&u);CHKERRQ(ierr);
207 ierr = PetscLogStagePop();CHKERRQ(ierr);
208 }
209 for (r = 1; r <= Nr; ++r) {
210 ierr = DMDestroy(&dm[r]);CHKERRQ(ierr);
211 }
212 /* Fit convergence rate */
213 ierr = PetscMalloc2(Nr+1, &x, Nr+1, &y);CHKERRQ(ierr);
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 ierr = PetscLinearRegression(Nr+1, x, y, &slope, &intercept);CHKERRQ(ierr);
220 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
221 alpha[f] = -slope * dim;
222 }
223 ierr = PetscFree2(x, y);CHKERRQ(ierr);
224 ierr = PetscFree2(dm, dof);CHKERRQ(ierr);
225 /* Restore solver */
226 ierr = TSReset(ts);CHKERRQ(ierr);
227 {
228 /* PCReset() does not wipe out the level structure */
229 SNES snes;
230 KSP ksp;
231 PC pc;
232
233 ierr = TSGetSNES(ts, &snes);CHKERRQ(ierr);
234 ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
235 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
236 ierr = PCMGSetLevels(pc, oldnlev, NULL);CHKERRQ(ierr);
237 ierr = DMSetRefineLevel(ce->idm, oldlevel);CHKERRQ(ierr); /* The damn DMCoarsen() calls in PCMG can reset this */
238 }
239 ierr = TSSetDM(ts, ce->idm);CHKERRQ(ierr);
240 ierr = DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx);CHKERRQ(ierr);
241 ierr = DMTSSetIFunctionLocal(ce->idm, DMPlexTSComputeIFunctionFEM, ctx);CHKERRQ(ierr);
242 ierr = DMTSSetIJacobianLocal(ce->idm, DMPlexTSComputeIJacobianFEM, ctx);CHKERRQ(ierr);
243 ierr = TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);CHKERRQ(ierr);
244 ierr = TSSetTime(ts, 0.0);CHKERRQ(ierr);
245 ierr = TSSetStepNumber(ts, 0);CHKERRQ(ierr);
246 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
247 ierr = TSSetSolution(ts, uInitial);CHKERRQ(ierr);
248 ierr = PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial);CHKERRQ(ierr);
249 PetscFunctionReturn(0);
250 }
251
PetscConvEstUseTS(PetscConvEst ce,PetscBool checkTemporal)252 PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
253 {
254 PetscFunctionBegin;
255 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1);
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 PetscFunctionReturn(0);
265 }
266