xref: /petsc/src/ts/utils/tsconvest.c (revision bef158480efac06de457f7a665168877ab3c2fd7) 
1 #include <petscconvest.h>            /*I "petscconvest.h" I*/
2 #include <petscts.h>
3 #include <petscdmplex.h>
4 
5 #include <petsc/private/petscconvestimpl.h>
6 
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 
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 
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 
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 
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 
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