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: +  sname - name of user-defined adaptivity scheme
 32: -  function - routine to create method context

 34:    Level: advanced

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

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

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

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

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

 62:   Not Collective

 64:   Level: advanced

 66: .seealso: [](ch_ts), `TSGLLEAdapt`, `TSGLLE`, `TSGLLEAdaptRegisterDestroy()`
 67: @*/
 68: PetscErrorCode TSGLLEAdaptRegisterAll(void)
 69: {
 70:   PetscFunctionBegin;
 71:   if (TSGLLEAdaptRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
 72:   TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
 73:   PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_NONE, TSGLLEAdaptCreate_None));
 74:   PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_SIZE, TSGLLEAdaptCreate_Size));
 75:   PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_BOTH, TSGLLEAdaptCreate_Both));
 76:   PetscFunctionReturn(PETSC_SUCCESS);
 77: }

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

 83:   Level: developer

 85: .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEAdapt`, `TSGLLEAdaptInitializePackage()`
 86: @*/
 87: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
 88: {
 89:   PetscFunctionBegin;
 90:   PetscCall(PetscFunctionListDestroy(&TSGLLEAdaptList));
 91:   TSGLLEAdaptPackageInitialized = PETSC_FALSE;
 92:   TSGLLEAdaptRegisterAllCalled  = PETSC_FALSE;
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

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

100:   Level: developer

102: .seealso: [](ch_ts), `PetscInitialize()`, `TSGLLEAdapt`, `TSGLLEAdaptFinalizePackage()`
103: @*/
104: PetscErrorCode TSGLLEAdaptInitializePackage(void)
105: {
106:   PetscFunctionBegin;
107:   if (TSGLLEAdaptPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
108:   TSGLLEAdaptPackageInitialized = PETSC_TRUE;
109:   PetscCall(PetscClassIdRegister("TSGLLEAdapt", &TSGLLEADAPT_CLASSID));
110:   PetscCall(TSGLLEAdaptRegisterAll());
111:   PetscCall(PetscRegisterFinalize(TSGLLEAdaptFinalizePackage));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
116: {
117:   PetscErrorCode (*r)(TSGLLEAdapt);

119:   PetscFunctionBegin;
120:   PetscCall(PetscFunctionListFind(TSGLLEAdaptList, type, &r));
121:   PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAdapt type \"%s\" given", type);
122:   if (((PetscObject)adapt)->type_name) PetscUseTypeMethod(adapt, destroy);
123:   PetscCall((*r)(adapt));
124:   PetscCall(PetscObjectChangeTypeName((PetscObject)adapt, type));
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
129: {
130:   PetscFunctionBegin;
131:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix));
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }

135: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
136: {
137:   PetscBool iascii;

139:   PetscFunctionBegin;
140:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
141:   if (iascii) {
142:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adapt, viewer));
143:     if (adapt->ops->view) {
144:       PetscCall(PetscViewerASCIIPushTab(viewer));
145:       PetscUseTypeMethod(adapt, view, viewer);
146:       PetscCall(PetscViewerASCIIPopTab(viewer));
147:     }
148:   }
149:   PetscFunctionReturn(PETSC_SUCCESS);
150: }

152: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
153: {
154:   PetscFunctionBegin;
155:   if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS);
157:   if (--((PetscObject)(*adapt))->refct > 0) {
158:     *adapt = NULL;
159:     PetscFunctionReturn(PETSC_SUCCESS);
160:   }
161:   PetscTryTypeMethod((*adapt), destroy);
162:   PetscCall(PetscHeaderDestroy(adapt));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems *PetscOptionsObject)
167: {
168:   char      type[256] = TSGLLEADAPT_BOTH;
169:   PetscBool flg;

171:   PetscFunctionBegin;
172:   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
173:   * function can only be called from inside TSSetFromOptions_GLLE()  */
174:   PetscOptionsHeadBegin(PetscOptionsObject, "TSGLLE Adaptivity options");
175:   PetscCall(PetscOptionsFList("-ts_adapt_type", "Algorithm to use for adaptivity", "TSGLLEAdaptSetType", TSGLLEAdaptList, ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type, type, sizeof(type), &flg));
176:   if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSGLLEAdaptSetType(adapt, type));
177:   PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject);
178:   PetscOptionsHeadEnd();
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }

182: 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)
183: {
184:   PetscFunctionBegin;
192:   PetscUseTypeMethod(adapt, choose, n, orders, errors, cost, cur, h, tleft, next_sc, next_h, finish);
193:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

196: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
197: {
198:   TSGLLEAdapt adapt;

200:   PetscFunctionBegin;
201:   *inadapt = NULL;
202:   PetscCall(PetscHeaderCreate(adapt, TSGLLEADAPT_CLASSID, "TSGLLEAdapt", "General Linear adaptivity", "TS", comm, TSGLLEAdaptDestroy, TSGLLEAdaptView));
203:   *inadapt = adapt;
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

207: /*
208:    Implementations
209: */

211: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
212: {
213:   PetscFunctionBegin;
214:   PetscCall(PetscFree(adapt->data));
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: /* -------------------------------- None ----------------------------------- */
219: typedef struct {
220:   PetscInt  scheme;
221:   PetscReal h;
222: } TSGLLEAdapt_None;

224: 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)
225: {
226:   PetscFunctionBegin;
227:   *next_sc = cur;
228:   *next_h  = h;
229:   if (*next_h > tleft) {
230:     *finish = PETSC_TRUE;
231:     *next_h = tleft;
232:   } else *finish = PETSC_FALSE;
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
237: {
238:   TSGLLEAdapt_None *a;

240:   PetscFunctionBegin;
241:   PetscCall(PetscNew(&a));
242:   adapt->data         = (void *)a;
243:   adapt->ops->choose  = TSGLLEAdaptChoose_None;
244:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: /* -------------------------------- Size ----------------------------------- */
249: typedef struct {
250:   PetscReal desired_h;
251: } TSGLLEAdapt_Size;

253: 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)
254: {
255:   TSGLLEAdapt_Size *sz  = (TSGLLEAdapt_Size *)adapt->data;
256:   PetscReal         dec = 0.2, inc = 5.0, safe = 0.9, optimal, last_desired_h;

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

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

270:   if (*next_h > tleft) {
271:     *finish = PETSC_TRUE;
272:     *next_h = tleft;
273:   } else *finish = PETSC_FALSE;
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }

277: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
278: {
279:   TSGLLEAdapt_Size *a;

281:   PetscFunctionBegin;
282:   PetscCall(PetscNew(&a));
283:   adapt->data         = (void *)a;
284:   adapt->ops->choose  = TSGLLEAdaptChoose_Size;
285:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

289: /* -------------------------------- Both ----------------------------------- */
290: typedef struct {
291:   PetscInt  count_at_order;
292:   PetscReal desired_h;
293: } TSGLLEAdapt_Both;

295: 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)
296: {
297:   TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both *)adapt->data;
298:   PetscReal         dec = 0.2, inc = 5.0, safe = 0.9;
299:   struct {
300:     PetscInt  id;
301:     PetscReal h, eff;
302:   } best = {-1, 0, 0}, trial = {-1, 0, 0}, current = {-1, 0, 0};
303:   PetscInt i;

305:   PetscFunctionBegin;
306:   for (i = 0; i < n; i++) {
307:     PetscReal optimal;
308:     trial.id  = i;
309:     optimal   = PetscPowReal((PetscReal)errors[i], (PetscReal)-1. / (safe * orders[i]));
310:     trial.h   = h * optimal;
311:     trial.eff = trial.h / cost[i];
312:     if (trial.eff > best.eff) PetscCall(PetscArraycpy(&best, &trial, 1));
313:     if (i == cur) PetscCall(PetscArraycpy(&current, &trial, 1));
314:   }
315:   /* Only switch orders if the scheme offers significant benefits over the current one.
316:   When the scheme is not changing, only change step size if it offers significant benefits. */
317:   if (best.eff < 1.2 * current.eff || both->count_at_order < orders[cur] + 2) {
318:     PetscReal last_desired_h;
319:     *next_sc        = current.id;
320:     last_desired_h  = both->desired_h;
321:     both->desired_h = PetscMax(h * dec, PetscMin(h * inc, current.h));
322:     *next_h         = (both->count_at_order > 0) ? PetscSqrtReal(last_desired_h * both->desired_h) : both->desired_h;
323:     both->count_at_order++;
324:   } else {
325:     PetscReal rat        = cost[best.id] / cost[cur];
326:     *next_sc             = best.id;
327:     *next_h              = PetscMax(h * rat * dec, PetscMin(h * rat * inc, best.h));
328:     both->count_at_order = 0;
329:     both->desired_h      = best.h;
330:   }

332:   if (*next_h > tleft) {
333:     *finish = PETSC_TRUE;
334:     *next_h = tleft;
335:   } else *finish = PETSC_FALSE;
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
340: {
341:   TSGLLEAdapt_Both *a;

343:   PetscFunctionBegin;
344:   PetscCall(PetscNew(&a));
345:   adapt->data         = (void *)a;
346:   adapt->ops->choose  = TSGLLEAdaptChoose_Both;
347:   adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }