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:   PetscFunctionBegin;
 12:   PetscCall(PetscObjectGetClassId(ce->solver, &id));
 13:   PetscCheck(id == TS_CLASSID, PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS");
 14:   PetscCall(TSGetDM((TS)ce->solver, &ce->idm));
 15:   PetscFunctionReturn(PETSC_SUCCESS);
 16: }

 18: static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
 19: {
 20:   PetscFunctionBegin;
 21:   PetscCall(TSComputeInitialCondition((TS)ce->solver, u));
 22:   PetscFunctionReturn(PETSC_SUCCESS);
 23: }

 25: static PetscErrorCode PetscConvEstComputeErrorTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[])
 26: {
 27:   TS ts = (TS)ce->solver;
 28:   PetscErrorCode (*exactError)(TS, Vec, Vec);

 30:   PetscFunctionBegin;
 31:   PetscCall(TSGetComputeExactError(ts, &exactError));
 32:   if (exactError) {
 33:     Vec      e;
 34:     PetscInt f;

 36:     PetscCall(VecDuplicate(u, &e));
 37:     PetscCall(TSComputeExactError(ts, u, e));
 38:     PetscCall(VecNorm(e, NORM_2, errors));
 39:     for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
 40:     PetscCall(VecDestroy(&e));
 41:   } else {
 42:     PetscReal t;

 44:     PetscCall(TSGetSolveTime(ts, &t));
 45:     PetscCall(DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors));
 46:   }
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: static PetscErrorCode PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce, PetscReal alpha[])
 51: {
 52:   TS         ts = (TS)ce->solver;
 53:   Vec        u, u0;
 54:   PetscReal *dt, *x, *y, slope, intercept;
 55:   PetscInt   Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r;

 57:   PetscFunctionBegin;
 58:   PetscCall(PetscMalloc1(Nr + 1, &dt));
 59:   PetscCall(TSGetTimeStep(ts, &dt[0]));
 60:   PetscCall(TSGetMaxSteps(ts, &oNs));
 61:   PetscCall(TSGetSolution(ts, &u0));
 62:   PetscCall(PetscObjectReference((PetscObject)u0));
 63:   Ns = oNs;
 64:   for (r = 0; r <= Nr; ++r) {
 65:     if (r > 0) {
 66:       dt[r] = dt[r - 1] / ce->r;
 67:       Ns    = PetscCeilReal(Ns * ce->r);
 68:     }
 69:     PetscCall(TSSetTime(ts, 0.0));
 70:     PetscCall(TSSetStepNumber(ts, 0));
 71:     PetscCall(TSSetTimeStep(ts, dt[r]));
 72:     PetscCall(TSSetMaxSteps(ts, Ns));
 73:     PetscCall(TSGetSolution(ts, &u));
 74:     PetscCall(PetscConvEstComputeInitialGuess(ce, r, NULL, u));
 75:     PetscCall(TSSolve(ts, NULL));
 76:     PetscCall(TSGetSolution(ts, &u));
 77:     PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
 78:     PetscCall(PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r * Nf]));
 79:     PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
 80:     for (f = 0; f < Nf; ++f) {
 81:       ce->dofs[r * Nf + f] = 1.0 / dt[r];
 82:       PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]));
 83:       PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]));
 84:     }
 85:     /* Monitor */
 86:     PetscCall(PetscConvEstMonitorDefault(ce, r));
 87:   }
 88:   /* Fit convergence rate */
 89:   if (Nr) {
 90:     PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
 91:     for (f = 0; f < Nf; ++f) {
 92:       for (r = 0; r <= Nr; ++r) {
 93:         x[r] = PetscLog10Real(dt[r]);
 94:         y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
 95:       }
 96:       PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
 97:       /* Since lg err = s lg dt + b */
 98:       alpha[f] = slope;
 99:     }
100:     PetscCall(PetscFree2(x, y));
101:   }
102:   /* Reset solver */
103:   PetscCall(TSReset(ts));
104:   PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
105:   PetscCall(TSSetTime(ts, 0.0));
106:   PetscCall(TSSetStepNumber(ts, 0));
107:   PetscCall(TSSetTimeStep(ts, dt[0]));
108:   PetscCall(TSSetMaxSteps(ts, oNs));
109:   PetscCall(TSSetSolution(ts, u0));
110:   PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, u0));
111:   PetscCall(VecDestroy(&u0));
112:   PetscCall(PetscFree(dt));
113:   PetscFunctionReturn(PETSC_SUCCESS);
114: }

116: static PetscErrorCode PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce, PetscReal alpha[])
117: {
118:   TS          ts = (TS)ce->solver;
119:   Vec         uInitial;
120:   DM         *dm;
121:   PetscObject disc;
122:   PetscReal  *x, *y, slope, intercept;
123:   PetscInt    Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev;
124:   PetscErrorCode (*ifunc)(DM, PetscReal, Vec, Vec, Vec, void *);
125:   PetscErrorCode (*ijac)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
126:   PetscErrorCode (*rhsfunc)(DM, PetscReal, Vec, Vec, void *);
127:   void *fctx, *jctx, *rctx;
128:   void *ctx;

130:   PetscFunctionBegin;
131:   PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r);
132:   PetscCall(DMGetDimension(ce->idm, &dim));
133:   PetscCall(DMGetApplicationContext(ce->idm, &ctx));
134:   PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE));
135:   PetscCall(DMGetRefineLevel(ce->idm, &oldlevel));
136:   PetscCall(PetscMalloc1((Nr + 1), &dm));
137:   PetscCall(TSGetSolution(ts, &uInitial));
138:   /* Loop over meshes */
139:   dm[0] = ce->idm;
140:   for (r = 0; r <= Nr; ++r) {
141:     Vec u;
142: #if defined(PETSC_USE_LOG)
143:     PetscLogStage stage;
144: #endif
145:     char        stageName[PETSC_MAX_PATH_LEN];
146:     const char *dmname, *uname;

148:     PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r));
149: #if defined(PETSC_USE_LOG)
150:     PetscCall(PetscLogStageGetId(stageName, &stage));
151:     if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage));
152: #endif
153:     PetscCall(PetscLogStagePush(stage));
154:     if (r > 0) {
155:       if (!ce->noRefine) {
156:         PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r]));
157:         PetscCall(DMSetCoarseDM(dm[r], dm[r - 1]));
158:       } else {
159:         DM cdm, rcdm;

161:         PetscCall(DMClone(dm[r - 1], &dm[r]));
162:         PetscCall(DMCopyDisc(dm[r - 1], dm[r]));
163:         PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm));
164:         PetscCall(DMGetCoordinateDM(dm[r], &rcdm));
165:         PetscCall(DMCopyDisc(cdm, rcdm));
166:       }
167:       PetscCall(DMCopyTransform(ce->idm, dm[r]));
168:       PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname));
169:       PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname));
170:       for (f = 0; f <= Nf; ++f) {
171:         PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);

173:         PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr));
174:         PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr));
175:       }
176:     }
177:     PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view"));
178:     /* Create solution */
179:     PetscCall(DMCreateGlobalVector(dm[r], &u));
180:     PetscCall(DMGetField(dm[r], 0, NULL, &disc));
181:     PetscCall(PetscObjectGetName(disc, &uname));
182:     PetscCall(PetscObjectSetName((PetscObject)u, uname));
183:     /* Setup solver */
184:     PetscCall(TSReset(ts));
185:     PetscCall(TSSetDM(ts, dm[r]));
186:     PetscCall(DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx));
187:     if (r > 0) {
188:       PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx));
189:       PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx));
190:       PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx));
191:       if (ifunc) PetscCall(DMTSSetIFunctionLocal(dm[r], ifunc, fctx));
192:       if (ijac) PetscCall(DMTSSetIJacobianLocal(dm[r], ijac, jctx));
193:       if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(dm[r], rhsfunc, rctx));
194:     }
195:     PetscCall(TSSetTime(ts, 0.0));
196:     PetscCall(TSSetStepNumber(ts, 0));
197:     PetscCall(TSSetFromOptions(ts));
198:     PetscCall(TSSetSolution(ts, u));
199:     PetscCall(VecDestroy(&u));
200:     /* Create initial guess */
201:     PetscCall(TSGetSolution(ts, &u));
202:     PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u));
203:     PetscCall(TSSolve(ts, NULL));
204:     PetscCall(TSGetSolution(ts, &u));
205:     PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0));
206:     PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * Nf]));
207:     PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0));
208:     for (f = 0; f < Nf; ++f) {
209:       PetscSection s, fs;
210:       PetscInt     lsize;

212:       /* Could use DMGetOutputDM() to add in Dirichlet dofs */
213:       PetscCall(DMGetLocalSection(dm[r], &s));
214:       PetscCall(PetscSectionGetField(s, f, &fs));
215:       PetscCall(PetscSectionGetConstrainedStorageSize(fs, &lsize));
216:       PetscCall(MPIU_Allreduce(&lsize, &ce->dofs[r * Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)ts)));
217:       PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f]));
218:       PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f]));
219:     }
220:     /* Monitor */
221:     PetscCall(PetscConvEstMonitorDefault(ce, r));
222:     if (!r) {
223:       /* PCReset() does not wipe out the level structure */
224:       SNES snes;
225:       KSP  ksp;
226:       PC   pc;

228:       PetscCall(TSGetSNES(ts, &snes));
229:       PetscCall(SNESGetKSP(snes, &ksp));
230:       PetscCall(KSPGetPC(ksp, &pc));
231:       PetscCall(PCMGGetLevels(pc, &oldnlev));
232:     }
233:     /* Cleanup */
234:     PetscCall(PetscLogStagePop());
235:   }
236:   PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx));
237:   PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx));
238:   PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx));
239:   for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r]));
240:   /* Fit convergence rate */
241:   PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y));
242:   for (f = 0; f < Nf; ++f) {
243:     for (r = 0; r <= Nr; ++r) {
244:       x[r] = PetscLog10Real(ce->dofs[r * Nf + f]);
245:       y[r] = PetscLog10Real(ce->errors[r * Nf + f]);
246:     }
247:     PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept));
248:     /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
249:     alpha[f] = -slope * dim;
250:   }
251:   PetscCall(PetscFree2(x, y));
252:   PetscCall(PetscFree(dm));
253:   /* Restore solver */
254:   PetscCall(TSReset(ts));
255:   {
256:     /* PCReset() does not wipe out the level structure */
257:     SNES snes;
258:     KSP  ksp;
259:     PC   pc;

261:     PetscCall(TSGetSNES(ts, &snes));
262:     PetscCall(SNESGetKSP(snes, &ksp));
263:     PetscCall(KSPGetPC(ksp, &pc));
264:     PetscCall(PCMGSetLevels(pc, oldnlev, NULL));
265:     PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */
266:   }
267:   PetscCall(TSSetDM(ts, ce->idm));
268:   PetscCall(DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx));
269:   if (ifunc) PetscCall(DMTSSetIFunctionLocal(ce->idm, ifunc, fctx));
270:   if (ijac) PetscCall(DMTSSetIJacobianLocal(ce->idm, ijac, jctx));
271:   if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(ce->idm, rhsfunc, rctx));
272:   PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING));
273:   PetscCall(TSSetTime(ts, 0.0));
274:   PetscCall(TSSetStepNumber(ts, 0));
275:   PetscCall(TSSetFromOptions(ts));
276:   PetscCall(TSSetSolution(ts, uInitial));
277:   PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial));
278:   PetscFunctionReturn(PETSC_SUCCESS);
279: }

281: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
282: {
283:   PetscFunctionBegin;
285:   ce->ops->setsolver    = PetscConvEstSetTS_Private;
286:   ce->ops->initguess    = PetscConvEstInitGuessTS_Private;
287:   ce->ops->computeerror = PetscConvEstComputeErrorTS_Private;
288:   if (checkTemporal) {
289:     ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
290:   } else {
291:     ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
292:   }
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }