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