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;
51: PetscReal *dt, *x, *y, slope, intercept;
52: PetscInt Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r;
54: TSGetSolution(ts, &u);
55: PetscMalloc1(Nr+1, &dt);
56: TSGetTimeStep(ts, &dt[0]);
57: TSGetMaxSteps(ts, &oNs);
58: Ns = oNs;
59: for (r = 0; r <= Nr; ++r) {
60: if (r > 0) {
61: dt[r] = dt[r-1]/ce->r;
62: Ns = PetscCeilReal(Ns*ce->r);
63: }
64: TSSetTime(ts, 0.0);
65: TSSetStepNumber(ts, 0);
66: TSSetTimeStep(ts, dt[r]);
67: TSSetMaxSteps(ts, Ns);
68: PetscConvEstComputeInitialGuess(ce, r, NULL, u);
69: TSSolve(ts, u);
70: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
71: PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r*Nf]);
72: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
73: for (f = 0; f < Nf; ++f) {
74: ce->dofs[r*Nf+f] = 1.0/dt[r];
75: PetscLogEventSetDof(ce->event, f, ce->dofs[r*Nf+f]);
76: PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
77: }
78: /* Monitor */
79: PetscConvEstMonitorDefault(ce, r);
80: }
81: /* Fit convergence rate */
82: if (Nr) {
83: PetscMalloc2(Nr+1, &x, Nr+1, &y);
84: for (f = 0; f < Nf; ++f) {
85: for (r = 0; r <= Nr; ++r) {
86: x[r] = PetscLog10Real(dt[r]);
87: y[r] = PetscLog10Real(ce->errors[r*Nf+f]);
88: }
89: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
90: /* Since lg err = s lg dt + b */
91: alpha[f] = slope;
92: }
93: PetscFree2(x, y);
94: }
95: /* Reset solver */
96: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
97: TSSetTime(ts, 0.0);
98: TSSetStepNumber(ts, 0);
99: TSSetTimeStep(ts, dt[0]);
100: TSSetMaxSteps(ts, oNs);
101: PetscConvEstComputeInitialGuess(ce, 0, NULL, u);
102: PetscFree(dt);
103: return 0;
104: }
106: static PetscErrorCode PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce, PetscReal alpha[])
107: {
108: TS ts = (TS) ce->solver;
109: Vec uInitial;
110: DM *dm;
111: PetscObject disc;
112: PetscReal *x, *y, slope, intercept;
113: PetscInt Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev;
114: void *ctx;
117: DMGetDimension(ce->idm, &dim);
118: DMGetApplicationContext(ce->idm, &ctx);
119: DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
120: DMGetRefineLevel(ce->idm, &oldlevel);
121: PetscMalloc1((Nr+1), &dm);
122: TSGetSolution(ts, &uInitial);
123: /* Loop over meshes */
124: dm[0] = ce->idm;
125: for (r = 0; r <= Nr; ++r) {
126: Vec u;
127: #if defined(PETSC_USE_LOG)
128: PetscLogStage stage;
129: #endif
130: char stageName[PETSC_MAX_PATH_LEN];
131: const char *dmname, *uname;
133: PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
134: #if defined(PETSC_USE_LOG)
135: PetscLogStageGetId(stageName, &stage);
136: if (stage < 0) PetscLogStageRegister(stageName, &stage);
137: #endif
138: PetscLogStagePush(stage);
139: if (r > 0) {
140: if (!ce->noRefine) {
141: DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
142: DMSetCoarseDM(dm[r], dm[r-1]);
143: } else {
144: DM cdm, rcdm;
146: DMClone(dm[r-1], &dm[r]);
147: DMCopyDisc(dm[r-1], dm[r]);
148: DMGetCoordinateDM(dm[r-1], &cdm);
149: DMGetCoordinateDM(dm[r], &rcdm);
150: DMCopyDisc(cdm, rcdm);
151: }
152: DMCopyTransform(ce->idm, dm[r]);
153: PetscObjectGetName((PetscObject) dm[r-1], &dmname);
154: PetscObjectSetName((PetscObject) dm[r], dmname);
155: for (f = 0; f <= Nf; ++f) {
156: PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);
158: DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
159: DMSetNullSpaceConstructor(dm[r], f, nspconstr);
160: }
161: }
162: DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
163: /* Create solution */
164: DMCreateGlobalVector(dm[r], &u);
165: DMGetField(dm[r], 0, NULL, &disc);
166: PetscObjectGetName(disc, &uname);
167: PetscObjectSetName((PetscObject) u, uname);
168: /* Setup solver */
169: TSReset(ts);
170: TSSetDM(ts, dm[r]);
171: DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx);
172: DMTSSetIFunctionLocal(dm[r], DMPlexTSComputeIFunctionFEM, ctx);
173: DMTSSetIJacobianLocal(dm[r], DMPlexTSComputeIJacobianFEM, ctx);
174: TSSetTime(ts, 0.0);
175: TSSetStepNumber(ts, 0);
176: TSSetFromOptions(ts);
177: /* Create initial guess */
178: PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
179: TSSolve(ts, u);
180: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
181: PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*Nf]);
182: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
183: for (f = 0; f < Nf; ++f) {
184: PetscSection s, fs;
185: PetscInt lsize;
187: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
188: DMGetLocalSection(dm[r], &s);
189: PetscSectionGetField(s, f, &fs);
190: PetscSectionGetConstrainedStorageSize(fs, &lsize);
191: MPI_Allreduce(&lsize, &ce->dofs[r*Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ts));
192: PetscLogEventSetDof(ce->event, f, ce->dofs[r*Nf+f]);
193: PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
194: }
195: /* Monitor */
196: PetscConvEstMonitorDefault(ce, r);
197: if (!r) {
198: /* PCReset() does not wipe out the level structure */
199: SNES snes;
200: KSP ksp;
201: PC pc;
203: TSGetSNES(ts, &snes);
204: SNESGetKSP(snes, &ksp);
205: KSPGetPC(ksp, &pc);
206: PCMGGetLevels(pc, &oldnlev);
207: }
208: /* Cleanup */
209: VecDestroy(&u);
210: PetscLogStagePop();
211: }
212: for (r = 1; r <= Nr; ++r) {
213: DMDestroy(&dm[r]);
214: }
215: /* Fit convergence rate */
216: PetscMalloc2(Nr+1, &x, Nr+1, &y);
217: for (f = 0; f < Nf; ++f) {
218: for (r = 0; r <= Nr; ++r) {
219: x[r] = PetscLog10Real(ce->dofs[r*Nf+f]);
220: y[r] = PetscLog10Real(ce->errors[r*Nf+f]);
221: }
222: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
223: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
224: alpha[f] = -slope * dim;
225: }
226: PetscFree2(x, y);
227: PetscFree(dm);
228: /* Restore solver */
229: TSReset(ts);
230: {
231: /* PCReset() does not wipe out the level structure */
232: SNES snes;
233: KSP ksp;
234: PC pc;
236: TSGetSNES(ts, &snes);
237: SNESGetKSP(snes, &ksp);
238: KSPGetPC(ksp, &pc);
239: PCMGSetLevels(pc, oldnlev, NULL);
240: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
241: }
242: TSSetDM(ts, ce->idm);
243: DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx);
244: DMTSSetIFunctionLocal(ce->idm, DMPlexTSComputeIFunctionFEM, ctx);
245: DMTSSetIJacobianLocal(ce->idm, DMPlexTSComputeIJacobianFEM, ctx);
246: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
247: TSSetTime(ts, 0.0);
248: TSSetStepNumber(ts, 0);
249: TSSetFromOptions(ts);
250: TSSetSolution(ts, uInitial);
251: PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial);
252: return 0;
253: }
255: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
256: {
258: ce->ops->setsolver = PetscConvEstSetTS_Private;
259: ce->ops->initguess = PetscConvEstInitGuessTS_Private;
260: ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
261: if (checkTemporal) {
262: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
263: } else {
264: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
265: }
266: return 0;
267: }