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
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: $ TSGLLEAdaptSetType(ts, "my_scheme")
45: or at runtime via the option
46: $ -ts_adapt_type my_scheme
48: .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegisterAll()`
49: @*/
50: PetscErrorCode TSGLLEAdaptRegister(const char sname[], PetscErrorCode (*function)(TSGLLEAdapt))
51: {
52: PetscFunctionBegin;
53: PetscCall(TSGLLEAdaptInitializePackage());
54: PetscCall(PetscFunctionListAdd(&TSGLLEAdaptList, sname, function));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
58: /*@C
59: TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in `TSGLLEAdapt`
61: Not Collective
63: Level: advanced
65: .seealso: [](ch_ts), `TSGLLEAdapt`, `TSGLLE`, `TSGLLEAdaptRegisterDestroy()`
66: @*/
67: PetscErrorCode TSGLLEAdaptRegisterAll(void)
68: {
69: PetscFunctionBegin;
70: if (TSGLLEAdaptRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
71: TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
72: PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_NONE, TSGLLEAdaptCreate_None));
73: PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_SIZE, TSGLLEAdaptCreate_Size));
74: PetscCall(TSGLLEAdaptRegister(TSGLLEADAPT_BOTH, TSGLLEAdaptCreate_Both));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: /*@C
79: TSGLLEAdaptFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
80: called from `PetscFinalize()`.
82: Level: developer
84: .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEAdapt`, `TSGLLEAdaptInitializePackage()`
85: @*/
86: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
87: {
88: PetscFunctionBegin;
89: PetscCall(PetscFunctionListDestroy(&TSGLLEAdaptList));
90: TSGLLEAdaptPackageInitialized = PETSC_FALSE;
91: TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
95: /*@C
96: TSGLLEAdaptInitializePackage - This function initializes everything in the `TSGLLEAdapt` package. It is
97: called from `TSInitializePackage()`.
99: Level: developer
101: .seealso: [](ch_ts), `PetscInitialize()`, `TSGLLEAdapt`, `TSGLLEAdaptFinalizePackage()`
102: @*/
103: PetscErrorCode TSGLLEAdaptInitializePackage(void)
104: {
105: PetscFunctionBegin;
106: if (TSGLLEAdaptPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
107: TSGLLEAdaptPackageInitialized = PETSC_TRUE;
108: PetscCall(PetscClassIdRegister("TSGLLEAdapt", &TSGLLEADAPT_CLASSID));
109: PetscCall(TSGLLEAdaptRegisterAll());
110: PetscCall(PetscRegisterFinalize(TSGLLEAdaptFinalizePackage));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
115: {
116: PetscErrorCode (*r)(TSGLLEAdapt);
118: PetscFunctionBegin;
119: PetscCall(PetscFunctionListFind(TSGLLEAdaptList, type, &r));
120: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAdapt type \"%s\" given", type);
121: if (((PetscObject)adapt)->type_name) PetscUseTypeMethod(adapt, destroy);
122: PetscCall((*r)(adapt));
123: PetscCall(PetscObjectChangeTypeName((PetscObject)adapt, type));
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
128: {
129: PetscFunctionBegin;
130: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
135: {
136: PetscBool iascii;
138: PetscFunctionBegin;
139: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
140: if (iascii) {
141: PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adapt, viewer));
142: if (adapt->ops->view) {
143: PetscCall(PetscViewerASCIIPushTab(viewer));
144: PetscUseTypeMethod(adapt, view, viewer);
145: PetscCall(PetscViewerASCIIPopTab(viewer));
146: }
147: }
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
152: {
153: PetscFunctionBegin;
154: if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS);
156: if (--((PetscObject)*adapt)->refct > 0) {
157: *adapt = NULL;
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
160: PetscTryTypeMethod(*adapt, destroy);
161: PetscCall(PetscHeaderDestroy(adapt));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems *PetscOptionsObject)
166: {
167: char type[256] = TSGLLEADAPT_BOTH;
168: PetscBool flg;
170: PetscFunctionBegin;
171: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
172: * function can only be called from inside TSSetFromOptions_GLLE() */
173: PetscOptionsHeadBegin(PetscOptionsObject, "TSGLLE Adaptivity options");
174: 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));
175: if (flg || !((PetscObject)adapt)->type_name) PetscCall(TSGLLEAdaptSetType(adapt, type));
176: PetscTryTypeMethod(adapt, setfromoptions, PetscOptionsObject);
177: PetscOptionsHeadEnd();
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: 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)
182: {
183: PetscFunctionBegin;
185: PetscAssertPointer(orders, 3);
186: PetscAssertPointer(errors, 4);
187: PetscAssertPointer(cost, 5);
188: PetscAssertPointer(next_sc, 9);
189: PetscAssertPointer(next_h, 10);
190: PetscAssertPointer(finish, 11);
191: PetscUseTypeMethod(adapt, choose, n, orders, errors, cost, cur, h, tleft, next_sc, next_h, finish);
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
196: {
197: TSGLLEAdapt adapt;
199: PetscFunctionBegin;
200: *inadapt = NULL;
201: PetscCall(PetscHeaderCreate(adapt, TSGLLEADAPT_CLASSID, "TSGLLEAdapt", "General Linear adaptivity", "TS", comm, TSGLLEAdaptDestroy, TSGLLEAdaptView));
202: *inadapt = adapt;
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*
207: Implementations
208: */
210: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
211: {
212: PetscFunctionBegin;
213: PetscCall(PetscFree(adapt->data));
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /* -------------------------------- None ----------------------------------- */
218: typedef struct {
219: PetscInt scheme;
220: PetscReal h;
221: } TSGLLEAdapt_None;
223: 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)
224: {
225: PetscFunctionBegin;
226: *next_sc = cur;
227: *next_h = h;
228: if (*next_h > tleft) {
229: *finish = PETSC_TRUE;
230: *next_h = tleft;
231: } else *finish = PETSC_FALSE;
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
236: {
237: TSGLLEAdapt_None *a;
239: PetscFunctionBegin;
240: PetscCall(PetscNew(&a));
241: adapt->data = (void *)a;
242: adapt->ops->choose = TSGLLEAdaptChoose_None;
243: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /* -------------------------------- Size ----------------------------------- */
248: typedef struct {
249: PetscReal desired_h;
250: } TSGLLEAdapt_Size;
252: 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)
253: {
254: TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size *)adapt->data;
255: PetscReal dec = 0.2, inc = 5.0, safe = 0.9, optimal, last_desired_h;
257: PetscFunctionBegin;
258: *next_sc = cur;
259: optimal = PetscPowReal((PetscReal)errors[cur], (PetscReal)-1. / (safe * orders[cur]));
260: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
261: * one that would have been taken (without smoothing) on the last step. */
262: last_desired_h = sz->desired_h;
263: sz->desired_h = h * PetscMax(dec, PetscMin(inc, optimal)); /* Trim to [dec,inc] */
265: /* Normally only happens on the first step */
266: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
267: else *next_h = sz->desired_h;
269: if (*next_h > tleft) {
270: *finish = PETSC_TRUE;
271: *next_h = tleft;
272: } else *finish = PETSC_FALSE;
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
277: {
278: TSGLLEAdapt_Size *a;
280: PetscFunctionBegin;
281: PetscCall(PetscNew(&a));
282: adapt->data = (void *)a;
283: adapt->ops->choose = TSGLLEAdaptChoose_Size;
284: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /* -------------------------------- Both ----------------------------------- */
289: typedef struct {
290: PetscInt count_at_order;
291: PetscReal desired_h;
292: } TSGLLEAdapt_Both;
294: 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)
295: {
296: TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both *)adapt->data;
297: PetscReal dec = 0.2, inc = 5.0, safe = 0.9;
298: struct {
299: PetscInt id;
300: PetscReal h, eff;
301: } best = {-1, 0, 0}, trial = {-1, 0, 0}, current = {-1, 0, 0};
302: PetscInt i;
304: PetscFunctionBegin;
305: for (i = 0; i < n; i++) {
306: PetscReal optimal;
307: trial.id = i;
308: optimal = PetscPowReal((PetscReal)errors[i], (PetscReal)-1. / (safe * orders[i]));
309: trial.h = h * optimal;
310: trial.eff = trial.h / cost[i];
311: if (trial.eff > best.eff) PetscCall(PetscArraycpy(&best, &trial, 1));
312: if (i == cur) PetscCall(PetscArraycpy(¤t, &trial, 1));
313: }
314: /* Only switch orders if the scheme offers significant benefits over the current one.
315: When the scheme is not changing, only change step size if it offers significant benefits. */
316: if (best.eff < 1.2 * current.eff || both->count_at_order < orders[cur] + 2) {
317: PetscReal last_desired_h;
318: *next_sc = current.id;
319: last_desired_h = both->desired_h;
320: both->desired_h = PetscMax(h * dec, PetscMin(h * inc, current.h));
321: *next_h = (both->count_at_order > 0) ? PetscSqrtReal(last_desired_h * both->desired_h) : both->desired_h;
322: both->count_at_order++;
323: } else {
324: PetscReal rat = cost[best.id] / cost[cur];
325: *next_sc = best.id;
326: *next_h = PetscMax(h * rat * dec, PetscMin(h * rat * inc, best.h));
327: both->count_at_order = 0;
328: both->desired_h = best.h;
329: }
331: if (*next_h > tleft) {
332: *finish = PETSC_TRUE;
333: *next_h = tleft;
334: } else *finish = PETSC_FALSE;
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
339: {
340: TSGLLEAdapt_Both *a;
342: PetscFunctionBegin;
343: PetscCall(PetscNew(&a));
344: adapt->data = (void *)a;
345: adapt->ops->choose = TSGLLEAdaptChoose_Both;
346: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }