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: PetscObjectGetClassId(ce->solver, &id);
13: TSGetDM((TS)ce->solver, &ce->idm);
14: return 0;
15: }
17: static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
18: {
19: TSComputeInitialCondition((TS)ce->solver, u);
20: return 0;
21: }
23: static PetscErrorCode PetscConvEstComputeErrorTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
24: {
25: TS ts = (TS)ce->solver;
26: PetscErrorCode (*exactError)(TS, Vec, Vec);
28: TSGetComputeExactError(ts, &exactError);
29: if (exactError) {
30: Vec e;
31: PetscInt f;
33: VecDuplicate(u, &e);
34: TSComputeExactError(ts, u, e);
35: VecNorm(e, NORM_2, errors);
36: for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
37: VecDestroy(&e);
38: } else {
39: PetscReal t;
41: TSGetSolveTime(ts, &t);
42: DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors);
43: }
44: return 0;
45: }
47: static PetscErrorCode PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce, PetscReal alpha[])
48: {
49: TS ts = (TS)ce->solver;
50: Vec u, u0;
51: PetscReal *dt, *x, *y, slope, intercept;
52: PetscInt Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r;
54: PetscMalloc1(Nr + 1, &dt);
55: TSGetTimeStep(ts, &dt[0]);
56: TSGetMaxSteps(ts, &oNs);
57: TSGetSolution(ts, &u0);
58: PetscObjectReference((PetscObject)u0);
59: Ns = oNs;
60: for (r = 0; r <= Nr; ++r) {
61: if (r > 0) {
62: dt[r] = dt[r - 1] / ce->r;
63: Ns = PetscCeilReal(Ns * ce->r);
64: }
65: TSSetTime(ts, 0.0);
66: TSSetStepNumber(ts, 0);
67: TSSetTimeStep(ts, dt[r]);
68: TSSetMaxSteps(ts, Ns);
69: TSGetSolution(ts, &u);
70: PetscConvEstComputeInitialGuess(ce, r, NULL, u);
71: TSSolve(ts, NULL);
72: TSGetSolution(ts, &u);
73: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
74: PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r * Nf]);
75: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
76: for (f = 0; f < Nf; ++f) {
77: ce->dofs[r * Nf + f] = 1.0 / dt[r];
78: PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]);
79: PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]);
80: }
81: /* Monitor */
82: PetscConvEstMonitorDefault(ce, r);
83: }
84: /* Fit convergence rate */
85: if (Nr) {
86: PetscMalloc2(Nr + 1, &x, Nr + 1, &y);
87: for (f = 0; f < Nf; ++f) {
88: for (r = 0; r <= Nr; ++r) {
89: x[r] = PetscLog10Real(dt[r]);
90: y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
91: }
92: PetscLinearRegression(Nr + 1, x, y, &slope, &intercept);
93: /* Since lg err = s lg dt + b */
94: alpha[f] = slope;
95: }
96: PetscFree2(x, y);
97: }
98: /* Reset solver */
99: TSReset(ts);
100: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
101: TSSetTime(ts, 0.0);
102: TSSetStepNumber(ts, 0);
103: TSSetTimeStep(ts, dt[0]);
104: TSSetMaxSteps(ts, oNs);
105: TSSetSolution(ts, u0);
106: PetscConvEstComputeInitialGuess(ce, 0, NULL, u0);
107: VecDestroy(&u0);
108: PetscFree(dt);
109: return 0;
110: }
112: static PetscErrorCode PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce, PetscReal alpha[])
113: {
114: TS ts = (TS)ce->solver;
115: Vec uInitial;
116: DM *dm;
117: PetscObject disc;
118: PetscReal *x, *y, slope, intercept;
119: PetscInt Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev;
120: PetscErrorCode (*ifunc)(DM, PetscReal, Vec, Vec, Vec, void *);
121: PetscErrorCode (*ijac)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
122: PetscErrorCode (*rhsfunc)(DM, PetscReal, Vec, Vec, void *);
123: void *fctx, *jctx, *rctx;
124: void *ctx;
127: DMGetDimension(ce->idm, &dim);
128: DMGetApplicationContext(ce->idm, &ctx);
129: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
130: DMGetRefineLevel(ce->idm, &oldlevel);
131: PetscMalloc1((Nr + 1), &dm);
132: TSGetSolution(ts, &uInitial);
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;
143: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r);
144: #if defined(PETSC_USE_LOG)
145: PetscLogStageGetId(stageName, &stage);
146: if (stage < 0) PetscLogStageRegister(stageName, &stage);
147: #endif
148: PetscLogStagePush(stage);
149: if (r > 0) {
150: if (!ce->noRefine) {
151: DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]);
152: DMSetCoarseDM(dm[r], dm[r - 1]);
153: } else {
154: DM cdm, rcdm;
156: DMClone(dm[r - 1], &dm[r]);
157: DMCopyDisc(dm[r - 1], dm[r]);
158: DMGetCoordinateDM(dm[r - 1], &cdm);
159: DMGetCoordinateDM(dm[r], &rcdm);
160: DMCopyDisc(cdm, rcdm);
161: }
162: DMCopyTransform(ce->idm, dm[r]);
163: PetscObjectGetName((PetscObject)dm[r - 1], &dmname);
164: PetscObjectSetName((PetscObject)dm[r], dmname);
165: for (f = 0; f <= Nf; ++f) {
166: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
168: DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr);
169: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
170: }
171: }
172: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
173: /* Create solution */
174: DMCreateGlobalVector(dm[r], &u);
175: DMGetField(dm[r], 0, NULL, &disc);
176: PetscObjectGetName(disc, &uname);
177: PetscObjectSetName((PetscObject)u, uname);
178: /* Setup solver */
179: TSReset(ts);
180: TSSetDM(ts, dm[r]);
181: DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx);
182: if (r > 0) {
183: DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx);
184: DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx);
185: DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx);
186: if (ifunc) DMTSSetIFunctionLocal(dm[r], ifunc, fctx);
187: if (ijac) DMTSSetIJacobianLocal(dm[r], ijac, jctx);
188: if (rhsfunc) DMTSSetRHSFunctionLocal(dm[r], rhsfunc, rctx);
189: }
190: TSSetTime(ts, 0.0);
191: TSSetStepNumber(ts, 0);
192: TSSetFromOptions(ts);
193: TSSetSolution(ts, u);
194: VecDestroy(&u);
195: /* Create initial guess */
196: TSGetSolution(ts, &u);
197: PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
198: TSSolve(ts, NULL);
199: TSGetSolution(ts, &u);
200: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
201: PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * Nf]);
202: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
203: for (f = 0; f < Nf; ++f) {
204: PetscSection s, fs;
205: PetscInt lsize;
207: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
208: DMGetLocalSection(dm[r], &s);
209: PetscSectionGetField(s, f, &fs);
210: PetscSectionGetConstrainedStorageSize(fs, &lsize);
211: MPI_Allreduce(&lsize, &ce->dofs[r * Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)ts));
212: PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]);
213: PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]);
214: }
215: /* Monitor */
216: PetscConvEstMonitorDefault(ce, r);
217: if (!r) {
218: /* PCReset() does not wipe out the level structure */
219: SNES snes;
220: KSP ksp;
221: PC pc;
223: TSGetSNES(ts, &snes);
224: SNESGetKSP(snes, &ksp);
225: KSPGetPC(ksp, &pc);
226: PCMGGetLevels(pc, &oldnlev);
227: }
228: /* Cleanup */
229: PetscLogStagePop();
230: }
231: DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx);
232: DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx);
233: DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx);
234: for (r = 1; r <= Nr; ++r) DMDestroy(&dm[r]);
235: /* Fit convergence rate */
236: PetscMalloc2(Nr + 1, &x, Nr + 1, &y);
237: for (f = 0; f < Nf; ++f) {
238: for (r = 0; r <= Nr; ++r) {
239: x[r] = PetscLog10Real(ce->dofs[r * Nf + f]);
240: y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
241: }
242: PetscLinearRegression(Nr + 1, x, y, &slope, &intercept);
243: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
244: alpha[f] = -slope * dim;
245: }
246: PetscFree2(x, y);
247: PetscFree(dm);
248: /* Restore solver */
249: TSReset(ts);
250: {
251: /* PCReset() does not wipe out the level structure */
252: SNES snes;
253: KSP ksp;
254: PC pc;
256: TSGetSNES(ts, &snes);
257: SNESGetKSP(snes, &ksp);
258: KSPGetPC(ksp, &pc);
259: PCMGSetLevels(pc, oldnlev, NULL);
260: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
261: }
262: TSSetDM(ts, ce->idm);
263: DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx);
264: if (ifunc) DMTSSetIFunctionLocal(ce->idm, ifunc, fctx);
265: if (ijac) DMTSSetIJacobianLocal(ce->idm, ijac, jctx);
266: if (rhsfunc) DMTSSetRHSFunctionLocal(ce->idm, rhsfunc, rctx);
267: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
268: TSSetTime(ts, 0.0);
269: TSSetStepNumber(ts, 0);
270: TSSetFromOptions(ts);
271: TSSetSolution(ts, uInitial);
272: PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial);
273: return 0;
274: }
276: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
277: {
279: ce->ops->setsolver = PetscConvEstSetTS_Private;
280: ce->ops->initguess = PetscConvEstInitGuessTS_Private;
281: ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
282: if (checkTemporal) {
283: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
284: } else {
285: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
286: }
287: return 0;
288: }