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