Actual source code: glleadapt.c

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

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

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

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

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

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

 27:   Not Collective, No Fortran Support

 29:   Input Parameters:
 30: + sname    - name of user-defined adaptivity scheme
 31: - function - routine to create method context

 33:   Level: advanced

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

 38:   Example Usage:
 39: .vb
 40:    TSGLLEAdaptRegister("my_scheme", MySchemeCreate);
 41: .ve

 43:   Then, your scheme can be chosen with the procedural interface via
 44: .vb
 45:   TSGLLEAdaptSetType(ts, "my_scheme")
 46: .ve
 47:   or at runtime via the option
 48: .vb
 49:   -ts_adapt_type my_scheme
 50: .ve

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

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

 65:   Not Collective

 67:   Level: advanced

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

 82: /*@C
 83:   TSGLLEAdaptFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
 84:   called from `PetscFinalize()`.

 86:   Level: developer

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

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

103:   Level: developer

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

118: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
119: {
120:   PetscErrorCode (*r)(TSGLLEAdapt);

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

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

138: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
139: {
140:   PetscBool iascii;

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

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

169: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems PetscOptionsObject)
170: {
171:   char      type[256] = TSGLLEADAPT_BOTH;
172:   PetscBool flg;

174:   PetscFunctionBegin;
175:   /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
176:   * function can only be called from inside TSSetFromOptions_GLLE()  */
177:   PetscOptionsHeadBegin(PetscOptionsObject, "TSGLLE Adaptivity options");
178:   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));
179:   if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSGLLEAdaptSetType(adapt, type));
180:   PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject);
181:   PetscOptionsHeadEnd();
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: 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)
186: {
187:   PetscFunctionBegin;
189:   PetscAssertPointer(orders, 3);
190:   PetscAssertPointer(errors, 4);
191:   PetscAssertPointer(cost, 5);
192:   PetscAssertPointer(next_sc, 9);
193:   PetscAssertPointer(next_h, 10);
194:   PetscAssertPointer(finish, 11);
195:   PetscUseTypeMethod(adapt, choose, n, orders, errors, cost, cur, h, tleft, next_sc, next_h, finish);
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
200: {
201:   TSGLLEAdapt adapt;

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

210: /*
211:    Implementations
212: */

214: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
215: {
216:   PetscFunctionBegin;
217:   PetscCall(PetscFree(adapt->data));
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: /* -------------------------------- None ----------------------------------- */
222: typedef struct {
223:   PetscInt  scheme;
224:   PetscReal h;
225: } TSGLLEAdapt_None;

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

239: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
240: {
241:   TSGLLEAdapt_None *a;

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

251: /* -------------------------------- Size ----------------------------------- */
252: typedef struct {
253:   PetscReal desired_h;
254: } TSGLLEAdapt_Size;

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

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

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

273:   if (*next_h > tleft) {
274:     *finish = PETSC_TRUE;
275:     *next_h = tleft;
276:   } else *finish = PETSC_FALSE;
277:   PetscFunctionReturn(PETSC_SUCCESS);
278: }

280: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
281: {
282:   TSGLLEAdapt_Size *a;

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

292: /* -------------------------------- Both ----------------------------------- */
293: typedef struct {
294:   PetscInt  count_at_order;
295:   PetscReal desired_h;
296: } TSGLLEAdapt_Both;

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

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

335:   if (*next_h > tleft) {
336:     *finish = PETSC_TRUE;
337:     *next_h = tleft;
338:   } else *finish = PETSC_FALSE;
339:   PetscFunctionReturn(PETSC_SUCCESS);
340: }

342: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
343: {
344:   TSGLLEAdapt_Both *a;

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