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(¤t, &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: }