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: /*@C
119: TSGLLEAdaptSetType - Sets the type of a `TSGLLEAdapt` step-size and order adaptivity object
121: Logically Collective
123: Input Parameters:
124: + adapt - the `TSGLLEAdapt` context
125: - type - the name of the adaptivity scheme, e.g. `TSGLLEADAPT_NONE`, `TSGLLEADAPT_SIZE`, `TSGLLEADAPT_BOTH`
127: Level: intermediate
129: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptCreate()`, `TSGLLEAdaptType`, `TSGLLEAdaptRegister()`
130: @*/
131: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
132: {
133: PetscErrorCode (*r)(TSGLLEAdapt);
135: PetscFunctionBegin;
136: PetscCall(PetscFunctionListFind(TSGLLEAdaptList, type, &r));
137: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAdapt type \"%s\" given", type);
138: if (((PetscObject)adapt)->type_name) PetscUseTypeMethod(adapt, destroy);
139: PetscCall((*r)(adapt));
140: PetscCall(PetscObjectChangeTypeName((PetscObject)adapt, type));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /*@
145: TSGLLEAdaptSetOptionsPrefix - Sets the prefix used for searching for `TSGLLEAdapt` options in the options database
147: Logically Collective
149: Input Parameters:
150: + adapt - the `TSGLLEAdapt` context
151: - prefix - the prefix to prepend to all option names
153: Level: advanced
155: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptSetFromOptions()`
156: @*/
157: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
158: {
159: PetscFunctionBegin;
160: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: /*@
165: TSGLLEAdaptView - Views a `TSGLLEAdapt` step-size and order adaptivity object
167: Collective
169: Input Parameters:
170: + adapt - the `TSGLLEAdapt` context
171: - viewer - the `PetscViewer` used to view the object
173: Level: intermediate
175: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptCreate()`, `PetscViewer`
176: @*/
177: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
178: {
179: PetscBool isascii;
181: PetscFunctionBegin;
182: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
183: if (isascii) {
184: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adapt, viewer));
185: if (adapt->ops->view) {
186: PetscCall(PetscViewerASCIIPushTab(viewer));
187: PetscUseTypeMethod(adapt, view, viewer);
188: PetscCall(PetscViewerASCIIPopTab(viewer));
189: }
190: }
191: PetscFunctionReturn(PETSC_SUCCESS);
192: }
194: /*@
195: TSGLLEAdaptDestroy - Destroys a `TSGLLEAdapt` context
197: Collective
199: Input Parameter:
200: . adapt - the `TSGLLEAdapt` context
202: Level: intermediate
204: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptCreate()`
205: @*/
206: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
207: {
208: PetscFunctionBegin;
209: if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS);
211: if (--((PetscObject)*adapt)->refct > 0) {
212: *adapt = NULL;
213: PetscFunctionReturn(PETSC_SUCCESS);
214: }
215: PetscTryTypeMethod(*adapt, destroy);
216: PetscCall(PetscHeaderDestroy(adapt));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@
221: TSGLLEAdaptSetFromOptions - Sets options from the options database for a `TSGLLEAdapt` context
223: Collective
225: Input Parameters:
226: + adapt - the `TSGLLEAdapt` context
227: - PetscOptionsObject - the `PetscOptionItems` used to process options
229: Options Database Key:
230: . -ts_adapt_type (none|size|both) - algorithm to use for adaptivity
232: Level: advanced
234: Note:
235: This function is currently intended for internal use from inside `TSSetFromOptions_GLLE()`.
237: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptSetType()`
238: @*/
239: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems PetscOptionsObject)
240: {
241: char type[256] = TSGLLEADAPT_BOTH;
242: PetscBool flg;
244: PetscFunctionBegin;
245: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
246: * function can only be called from inside TSSetFromOptions_GLLE() */
247: PetscOptionsHeadBegin(PetscOptionsObject, "TSGLLE Adaptivity options");
248: 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));
249: if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSGLLEAdaptSetType(adapt, type));
250: PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject);
251: PetscOptionsHeadEnd();
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /*@
256: TSGLLEAdaptChoose - Choose the next scheme and step size using a `TSGLLEAdapt` step-size and order controller
258: Logically Collective
260: Input Parameters:
261: + adapt - the `TSGLLEAdapt` context
262: . n - the number of candidate schemes
263: . orders - the orders of accuracy of the candidate schemes
264: . errors - the error estimates for each candidate scheme
265: . cost - the relative cost of each candidate scheme
266: . cur - the index of the currently active scheme
267: . h - the last step size that was taken
268: - tleft - the amount of remaining integration time
270: Output Parameters:
271: + next_sc - the index of the scheme to use next
272: . next_h - the step size to take next
273: - finish - `PETSC_TRUE` if `next_h` was truncated to `tleft` because the end of the interval has been reached
275: Level: developer
277: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptCreate()`, `TSGLLEAdaptSetType()`
278: @*/
279: 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)
280: {
281: PetscFunctionBegin;
283: PetscAssertPointer(orders, 3);
284: PetscAssertPointer(errors, 4);
285: PetscAssertPointer(cost, 5);
286: PetscAssertPointer(next_sc, 9);
287: PetscAssertPointer(next_h, 10);
288: PetscAssertPointer(finish, 11);
289: PetscUseTypeMethod(adapt, choose, n, orders, errors, cost, cur, h, tleft, next_sc, next_h, finish);
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: /*@
294: TSGLLEAdaptCreate - Create a `TSGLLEAdapt` step-size and order adaptivity object
296: Collective
298: Input Parameter:
299: . comm - the MPI communicator
301: Output Parameter:
302: . inadapt - the newly created `TSGLLEAdapt` context
304: Level: intermediate
306: Note:
307: Typically this is not called by users; use `TSGetAdapt()` on the enclosing `TS` instead.
309: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptSetType()`, `TSGLLEAdaptDestroy()`
310: @*/
311: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
312: {
313: TSGLLEAdapt adapt;
315: PetscFunctionBegin;
316: *inadapt = NULL;
317: PetscCall(PetscHeaderCreate(adapt, TSGLLEADAPT_CLASSID, "TSGLLEAdapt", "General Linear adaptivity", "TS", comm, TSGLLEAdaptDestroy, TSGLLEAdaptView));
318: *inadapt = adapt;
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
323: {
324: PetscFunctionBegin;
325: PetscCall(PetscFree(adapt->data));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: typedef struct {
330: PetscInt scheme;
331: PetscReal h;
332: } TSGLLEAdapt_None;
334: 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)
335: {
336: PetscFunctionBegin;
337: *next_sc = cur;
338: *next_h = h;
339: if (*next_h > tleft) {
340: *finish = PETSC_TRUE;
341: *next_h = tleft;
342: } else *finish = PETSC_FALSE;
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
347: {
348: TSGLLEAdapt_None *a;
350: PetscFunctionBegin;
351: PetscCall(PetscNew(&a));
352: adapt->data = (void *)a;
353: adapt->ops->choose = TSGLLEAdaptChoose_None;
354: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: typedef struct {
359: PetscReal desired_h;
360: } TSGLLEAdapt_Size;
362: 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)
363: {
364: TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size *)adapt->data;
365: PetscReal dec = 0.2, inc = 5.0, safe = 0.9, optimal, last_desired_h;
367: PetscFunctionBegin;
368: *next_sc = cur;
369: optimal = PetscPowReal((PetscReal)errors[cur], (PetscReal)-1. / (safe * orders[cur]));
370: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
371: * one that would have been taken (without smoothing) on the last step. */
372: last_desired_h = sz->desired_h;
373: sz->desired_h = h * PetscMax(dec, PetscMin(inc, optimal)); /* Trim to [dec,inc] */
375: /* Normally only happens on the first step */
376: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
377: else *next_h = sz->desired_h;
379: if (*next_h > tleft) {
380: *finish = PETSC_TRUE;
381: *next_h = tleft;
382: } else *finish = PETSC_FALSE;
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
387: {
388: TSGLLEAdapt_Size *a;
390: PetscFunctionBegin;
391: PetscCall(PetscNew(&a));
392: adapt->data = (void *)a;
393: adapt->ops->choose = TSGLLEAdaptChoose_Size;
394: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: typedef struct {
399: PetscInt count_at_order;
400: PetscReal desired_h;
401: } TSGLLEAdapt_Both;
403: 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)
404: {
405: TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both *)adapt->data;
406: PetscReal dec = 0.2, inc = 5.0, safe = 0.9;
407: struct {
408: PetscInt id;
409: PetscReal h, eff;
410: } best = {-1, 0, 0}, trial = {-1, 0, 0}, current = {-1, 0, 0};
411: PetscInt i;
413: PetscFunctionBegin;
414: for (i = 0; i < n; i++) {
415: PetscReal optimal;
416: trial.id = i;
417: optimal = PetscPowReal((PetscReal)errors[i], (PetscReal)-1. / (safe * orders[i]));
418: trial.h = h * optimal;
419: trial.eff = trial.h / cost[i];
420: if (trial.eff > best.eff) PetscCall(PetscArraycpy(&best, &trial, 1));
421: if (i == cur) PetscCall(PetscArraycpy(¤t, &trial, 1));
422: }
423: /* Only switch orders if the scheme offers significant benefits over the current one.
424: When the scheme is not changing, only change step size if it offers significant benefits. */
425: if (best.eff < 1.2 * current.eff || both->count_at_order < orders[cur] + 2) {
426: PetscReal last_desired_h;
427: *next_sc = current.id;
428: last_desired_h = both->desired_h;
429: both->desired_h = PetscMax(h * dec, PetscMin(h * inc, current.h));
430: *next_h = (both->count_at_order > 0) ? PetscSqrtReal(last_desired_h * both->desired_h) : both->desired_h;
431: both->count_at_order++;
432: } else {
433: PetscReal rat = cost[best.id] / cost[cur];
434: *next_sc = best.id;
435: *next_h = PetscMax(h * rat * dec, PetscMin(h * rat * inc, best.h));
436: both->count_at_order = 0;
437: both->desired_h = best.h;
438: }
440: if (*next_h > tleft) {
441: *finish = PETSC_TRUE;
442: *next_h = tleft;
443: } else *finish = PETSC_FALSE;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
448: {
449: TSGLLEAdapt_Both *a;
451: PetscFunctionBegin;
452: PetscCall(PetscNew(&a));
453: adapt->data = (void *)a;
454: adapt->ops->choose = TSGLLEAdaptChoose_Both;
455: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }