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;
13: PetscObjectGetClassId(ce->solver, &id);
14: if (id != TS_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS");
15: TSGetDM((TS) ce->solver, &ce->idm);
16: return(0);
17: }
19: static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
20: {
24: TSComputeInitialCondition((TS) ce->solver, u);
25: return(0);
26: }
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;
35: TSGetComputeExactError(ts, &exactError);
36: if (exactError) {
37: Vec e;
38: PetscInt f;
40: VecDuplicate(u, &e);
41: TSComputeExactError(ts, u, e);
42: VecNorm(e, NORM_2, errors);
43: for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
44: VecDestroy(&e);
45: } else {
46: PetscReal t;
48: TSGetSolveTime(ts, &t);
49: DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors);
50: }
51: return(0);
52: }
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;
63: TSGetSolution(ts, &u);
64: PetscMalloc1(Nr+1, &dt);
65: TSGetTimeStep(ts, &dt[0]);
66: TSGetMaxSteps(ts, &oNs);
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: TSSetTime(ts, 0.0);
74: TSSetStepNumber(ts, 0);
75: TSSetTimeStep(ts, dt[r]);
76: TSSetMaxSteps(ts, Ns);
77: PetscConvEstComputeInitialGuess(ce, r, NULL, u);
78: TSSolve(ts, u);
79: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
80: PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r*Nf]);
81: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
82: for (f = 0; f < Nf; ++f) {
83: PetscLogEventSetDof(ce->event, f, 1.0/dt[r]);
84: PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
85: }
86: /* Monitor */
87: PetscConvEstMonitorDefault(ce, r);
88: }
89: /* Fit convergence rate */
90: if (Nr) {
91: PetscMalloc2(Nr+1, &x, Nr+1, &y);
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: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
98: /* Since lg err = s lg dt + b */
99: alpha[f] = slope;
100: }
101: PetscFree2(x, y);
102: }
103: /* Reset solver */
104: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
105: TSSetTime(ts, 0.0);
106: TSSetStepNumber(ts, 0);
107: TSSetTimeStep(ts, dt[0]);
108: TSSetMaxSteps(ts, oNs);
109: PetscConvEstComputeInitialGuess(ce, 0, NULL, u);
110: PetscFree(dt);
111: return(0);
112: }
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 Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev;
122: void *ctx;
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: 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 %D", 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: DMTSSetIFunctionLocal(dm[r], DMPlexTSComputeIFunctionFEM, ctx);
183: DMTSSetIJacobianLocal(dm[r], DMPlexTSComputeIJacobianFEM, ctx);
184: TSSetTime(ts, 0.0);
185: TSSetStepNumber(ts, 0);
186: TSSetFromOptions(ts);
187: /* Create initial guess */
188: PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
189: TSSolve(ts, u);
190: PetscLogEventBegin(ce->event, ce, 0, 0, 0);
191: PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*Nf]);
192: PetscLogEventEnd(ce->event, ce, 0, 0, 0);
193: for (f = 0; f < Nf; ++f) {
194: PetscSection s, fs;
195: PetscInt lsize;
197: /* Could use DMGetOutputDM() to add in Dirichlet dofs */
198: DMGetLocalSection(dm[r], &s);
199: PetscSectionGetField(s, f, &fs);
200: PetscSectionGetConstrainedStorageSize(fs, &lsize);
201: MPI_Allreduce(&lsize, &ce->dofs[r*Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ts));
202: PetscLogEventSetDof(ce->event, f, ce->dofs[r*Nf+f]);
203: PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
204: }
205: /* Monitor */
206: PetscConvEstMonitorDefault(ce, r);
207: if (!r) {
208: /* PCReset() does not wipe out the level structure */
209: SNES snes;
210: KSP ksp;
211: PC pc;
213: TSGetSNES(ts, &snes);
214: SNESGetKSP(snes, &ksp);
215: KSPGetPC(ksp, &pc);
216: PCMGGetLevels(pc, &oldnlev);
217: }
218: /* Cleanup */
219: VecDestroy(&u);
220: PetscLogStagePop();
221: }
222: for (r = 1; r <= Nr; ++r) {
223: DMDestroy(&dm[r]);
224: }
225: /* Fit convergence rate */
226: PetscMalloc2(Nr+1, &x, Nr+1, &y);
227: for (f = 0; f < Nf; ++f) {
228: for (r = 0; r <= Nr; ++r) {
229: x[r] = PetscLog10Real(ce->dofs[r*Nf+f]);
230: y[r] = PetscLog10Real(ce->errors[r*Nf+f]);
231: }
232: PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
233: /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
234: alpha[f] = -slope * dim;
235: }
236: PetscFree2(x, y);
237: PetscFree(dm);
238: /* Restore solver */
239: TSReset(ts);
240: {
241: /* PCReset() does not wipe out the level structure */
242: SNES snes;
243: KSP ksp;
244: PC pc;
246: TSGetSNES(ts, &snes);
247: SNESGetKSP(snes, &ksp);
248: KSPGetPC(ksp, &pc);
249: PCMGSetLevels(pc, oldnlev, NULL);
250: DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
251: }
252: TSSetDM(ts, ce->idm);
253: DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx);
254: DMTSSetIFunctionLocal(ce->idm, DMPlexTSComputeIFunctionFEM, ctx);
255: DMTSSetIJacobianLocal(ce->idm, DMPlexTSComputeIJacobianFEM, ctx);
256: TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
257: TSSetTime(ts, 0.0);
258: TSSetStepNumber(ts, 0);
259: TSSetFromOptions(ts);
260: TSSetSolution(ts, uInitial);
261: PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial);
262: return(0);
263: }
265: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
266: {
269: ce->ops->setsolver = PetscConvEstSetTS_Private;
270: ce->ops->initguess = PetscConvEstInitGuessTS_Private;
271: ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
272: if (checkTemporal) {
273: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
274: } else {
275: ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
276: }
277: return(0);
278: }