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