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: }