Actual source code: glleadapt.c


  2: #include <../src/ts/impls/implicit/glle/glle.h>

  4: static PetscFunctionList TSGLLEAdaptList;
  5: static PetscBool         TSGLLEAdaptPackageInitialized;
  6: static PetscBool         TSGLLEAdaptRegisterAllCalled;
  7: static PetscClassId      TSGLLEADAPT_CLASSID;

  9: struct _TSGLLEAdaptOps {
 10:   PetscErrorCode (*choose)(TSGLLEAdapt, PetscInt, const PetscInt[], const PetscReal[], const PetscReal[], PetscInt, PetscReal, PetscReal, PetscInt *, PetscReal *, PetscBool *);
 11:   PetscErrorCode (*destroy)(TSGLLEAdapt);
 12:   PetscErrorCode (*view)(TSGLLEAdapt, PetscViewer);
 13:   PetscErrorCode (*setfromoptions)(TSGLLEAdapt, PetscOptionItems *);
 14: };

 16: struct _p_TSGLLEAdapt {
 17:   PETSCHEADER(struct _TSGLLEAdaptOps);
 18:   void *data;
 19: };

 21: PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt);
 22: PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt);
 23: PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt);

 25: /*@C
 26:    TSGLLEAdaptRegister -  adds a `TSGLLEAdapt` implementation

 28:    Not Collective

 30:    Input Parameters:
 31: +  name_scheme - name of user-defined adaptivity scheme
 32: -  routine_create - routine to create method context

 34:    Note:
 35:    `TSGLLEAdaptRegister()` may be called multiple times to add several user-defined families.

 37:    Sample usage:
 38: .vb
 39:    TSGLLEAdaptRegister("my_scheme",MySchemeCreate);
 40: .ve

 42:    Then, your scheme can be chosen with the procedural interface via
 43: $     TSGLLEAdaptSetType(ts,"my_scheme")
 44:    or at runtime via the option
 45: $     -ts_adapt_type my_scheme

 47:    Level: advanced

 49: .seealso: [](chapter_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegisterAll()`
 50: @*/
 51: PetscErrorCode TSGLLEAdaptRegister(const char sname[], PetscErrorCode (*function)(TSGLLEAdapt))
 52: {
 53:   TSGLLEAdaptInitializePackage();
 54:   PetscFunctionListAdd(&TSGLLEAdaptList, sname, function);
 55:   return 0;
 56: }

 58: /*@C
 59:   TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in `TSGLLEAdapt`

 61:   Not Collective

 63:   Level: advanced

 65: .seealso: [](chapter_ts), `TSGLLEAdapt`, `TSGLLE`, `TSGLLEAdaptRegisterDestroy()`
 66: @*/
 67: PetscErrorCode TSGLLEAdaptRegisterAll(void)
 68: {
 69:   if (TSGLLEAdaptRegisterAllCalled) return 0;
 70:   TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
 71:   TSGLLEAdaptRegister(TSGLLEADAPT_NONE, TSGLLEAdaptCreate_None);
 72:   TSGLLEAdaptRegister(TSGLLEADAPT_SIZE, TSGLLEAdaptCreate_Size);
 73:   TSGLLEAdaptRegister(TSGLLEADAPT_BOTH, TSGLLEAdaptCreate_Both);
 74:   return 0;
 75: }

 77: /*@C
 78:   TSGLLEFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
 79:   called from `PetscFinalize()`.

 81:   Level: developer

 83: .seealso: [](chapter_ts), `PetscFinalize()`, `TSGLLEAdapt`, `TSGLLEAdaptInitializePackage()`
 84: @*/
 85: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
 86: {
 87:   PetscFunctionListDestroy(&TSGLLEAdaptList);
 88:   TSGLLEAdaptPackageInitialized = PETSC_FALSE;
 89:   TSGLLEAdaptRegisterAllCalled  = PETSC_FALSE;
 90:   return 0;
 91: }

 93: /*@C
 94:   TSGLLEAdaptInitializePackage - This function initializes everything in the `TSGLLEAdapt` package. It is
 95:   called from `TSInitializePackage()`.

 97:   Level: developer

 99: .seealso: [](chapter_ts), `PetscInitialize()`, `TSGLLEAdapt`, `TSGLLEAdaptFinalizePackage()`
100: @*/
101: PetscErrorCode TSGLLEAdaptInitializePackage(void)
102: {
103:   if (TSGLLEAdaptPackageInitialized) return 0;
104:   TSGLLEAdaptPackageInitialized = PETSC_TRUE;
105:   PetscClassIdRegister("TSGLLEAdapt", &TSGLLEADAPT_CLASSID);
106:   TSGLLEAdaptRegisterAll();
107:   PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
108:   return 0;
109: }

111: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
112: {
113:   PetscErrorCode (*r)(TSGLLEAdapt);

115:   PetscFunctionListFind(TSGLLEAdaptList, type, &r);
117:   if (((PetscObject)adapt)->type_name) PetscUseTypeMethod(adapt, destroy);
118:   (*r)(adapt);
119:   PetscObjectChangeTypeName((PetscObject)adapt, type);
120:   return 0;
121: }

123: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
124: {
125:   PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix);
126:   return 0;
127: }

129: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
130: {
131:   PetscBool iascii;

133:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
134:   if (iascii) {
135:     PetscObjectPrintClassNamePrefixType((PetscObject)adapt, viewer);
136:     if (adapt->ops->view) {
137:       PetscViewerASCIIPushTab(viewer);
138:       PetscUseTypeMethod(adapt, view, viewer);
139:       PetscViewerASCIIPopTab(viewer);
140:     }
141:   }
142:   return 0;
143: }

145: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
146: {
147:   if (!*adapt) return 0;
149:   if (--((PetscObject)(*adapt))->refct > 0) {
150:     *adapt = NULL;
151:     return 0;
152:   }
153:   PetscTryTypeMethod((*adapt), destroy);
154:   PetscHeaderDestroy(adapt);
155:   return 0;
156: }

158: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems *PetscOptionsObject)
159: {
160:   char      type[256] = TSGLLEADAPT_BOTH;
161:   PetscBool flg;

163:   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
164:   * function can only be called from inside TSSetFromOptions_GLLE()  */
165:   PetscOptionsHeadBegin(PetscOptionsObject, "TSGLLE Adaptivity options");
166:   PetscOptionsFList("-ts_adapt_type", "Algorithm to use for adaptivity", "TSGLLEAdaptSetType", TSGLLEAdaptList, ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type, type, sizeof(type), &flg);
167:   if (flg || !((PetscObject)adapt)->type_name) TSGLLEAdaptSetType(adapt, type);
168:   PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject);
169:   PetscOptionsHeadEnd();
170:   return 0;
171: }

173: PetscErrorCode TSGLLEAdaptChoose(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish)
174: {
182:   PetscUseTypeMethod(adapt, choose, n, orders, errors, cost, cur, h, tleft, next_sc, next_h, finish);
183:   return 0;
184: }

186: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
187: {
188:   TSGLLEAdapt adapt;

190:   *inadapt = NULL;
191:   PetscHeaderCreate(adapt, TSGLLEADAPT_CLASSID, "TSGLLEAdapt", "General Linear adaptivity", "TS", comm, TSGLLEAdaptDestroy, TSGLLEAdaptView);
192:   *inadapt = adapt;
193:   return 0;
194: }

196: /*
197:    Implementations
198: */

200: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
201: {
202:   PetscFree(adapt->data);
203:   return 0;
204: }

206: /* -------------------------------- None ----------------------------------- */
207: typedef struct {
208:   PetscInt  scheme;
209:   PetscReal h;
210: } TSGLLEAdapt_None;

212: static PetscErrorCode TSGLLEAdaptChoose_None(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish)
213: {
214:   *next_sc = cur;
215:   *next_h  = h;
216:   if (*next_h > tleft) {
217:     *finish = PETSC_TRUE;
218:     *next_h = tleft;
219:   } else *finish = PETSC_FALSE;
220:   return 0;
221: }

223: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
224: {
225:   TSGLLEAdapt_None *a;

227:   PetscNew(&a);
228:   adapt->data         = (void *)a;
229:   adapt->ops->choose  = TSGLLEAdaptChoose_None;
230:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
231:   return 0;
232: }

234: /* -------------------------------- Size ----------------------------------- */
235: typedef struct {
236:   PetscReal desired_h;
237: } TSGLLEAdapt_Size;

239: static PetscErrorCode TSGLLEAdaptChoose_Size(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish)
240: {
241:   TSGLLEAdapt_Size *sz  = (TSGLLEAdapt_Size *)adapt->data;
242:   PetscReal         dec = 0.2, inc = 5.0, safe = 0.9, optimal, last_desired_h;

244:   *next_sc = cur;
245:   optimal  = PetscPowReal((PetscReal)errors[cur], (PetscReal)-1. / (safe * orders[cur]));
246:   /* Step sizes oscillate when there is no smoothing.  Here we use a geometric mean of the current step size and the
247:   * one that would have been taken (without smoothing) on the last step. */
248:   last_desired_h = sz->desired_h;
249:   sz->desired_h  = h * PetscMax(dec, PetscMin(inc, optimal)); /* Trim to [dec,inc] */

251:   /* Normally only happens on the first step */
252:   if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
253:   else *next_h = sz->desired_h;

255:   if (*next_h > tleft) {
256:     *finish = PETSC_TRUE;
257:     *next_h = tleft;
258:   } else *finish = PETSC_FALSE;
259:   return 0;
260: }

262: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
263: {
264:   TSGLLEAdapt_Size *a;

266:   PetscNew(&a);
267:   adapt->data         = (void *)a;
268:   adapt->ops->choose  = TSGLLEAdaptChoose_Size;
269:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
270:   return 0;
271: }

273: /* -------------------------------- Both ----------------------------------- */
274: typedef struct {
275:   PetscInt  count_at_order;
276:   PetscReal desired_h;
277: } TSGLLEAdapt_Both;

279: static PetscErrorCode TSGLLEAdaptChoose_Both(TSGLLEAdapt adapt, PetscInt n, const PetscInt orders[], const PetscReal errors[], const PetscReal cost[], PetscInt cur, PetscReal h, PetscReal tleft, PetscInt *next_sc, PetscReal *next_h, PetscBool *finish)
280: {
281:   TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both *)adapt->data;
282:   PetscReal         dec = 0.2, inc = 5.0, safe = 0.9;
283:   struct {
284:     PetscInt  id;
285:     PetscReal h, eff;
286:   } best = {-1, 0, 0}, trial = {-1, 0, 0}, current = {-1, 0, 0};
287:   PetscInt i;

289:   for (i = 0; i < n; i++) {
290:     PetscReal optimal;
291:     trial.id  = i;
292:     optimal   = PetscPowReal((PetscReal)errors[i], (PetscReal)-1. / (safe * orders[i]));
293:     trial.h   = h * optimal;
294:     trial.eff = trial.h / cost[i];
295:     if (trial.eff > best.eff) PetscArraycpy(&best, &trial, 1);
296:     if (i == cur) PetscArraycpy(&current, &trial, 1);
297:   }
298:   /* Only switch orders if the scheme offers significant benefits over the current one.
299:   When the scheme is not changing, only change step size if it offers significant benefits. */
300:   if (best.eff < 1.2 * current.eff || both->count_at_order < orders[cur] + 2) {
301:     PetscReal last_desired_h;
302:     *next_sc        = current.id;
303:     last_desired_h  = both->desired_h;
304:     both->desired_h = PetscMax(h * dec, PetscMin(h * inc, current.h));
305:     *next_h         = (both->count_at_order > 0) ? PetscSqrtReal(last_desired_h * both->desired_h) : both->desired_h;
306:     both->count_at_order++;
307:   } else {
308:     PetscReal rat        = cost[best.id] / cost[cur];
309:     *next_sc             = best.id;
310:     *next_h              = PetscMax(h * rat * dec, PetscMin(h * rat * inc, best.h));
311:     both->count_at_order = 0;
312:     both->desired_h      = best.h;
313:   }

315:   if (*next_h > tleft) {
316:     *finish = PETSC_TRUE;
317:     *next_h = tleft;
318:   } else *finish = PETSC_FALSE;
319:   return 0;
320: }

322: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
323: {
324:   TSGLLEAdapt_Both *a;

326:   PetscNew(&a);
327:   adapt->data         = (void *)a;
328:   adapt->ops->choose  = TSGLLEAdaptChoose_Both;
329:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
330:   return 0;
331: }