1 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/
2
3 static PetscFunctionList TSGLLEAdaptList;
4 static PetscBool TSGLLEAdaptPackageInitialized;
5 static PetscBool TSGLLEAdaptRegisterAllCalled;
6 static PetscClassId TSGLLEADAPT_CLASSID;
7
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 };
14
15 struct _p_TSGLLEAdapt {
16 PETSCHEADER(struct _TSGLLEAdaptOps);
17 void *data;
18 };
19
20 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt);
21 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt);
22 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt);
23
24 /*@C
25 TSGLLEAdaptRegister - adds a `TSGLLEAdapt` implementation
26
27 Not Collective, No Fortran Support
28
29 Input Parameters:
30 + sname - name of user-defined adaptivity scheme
31 - function - routine to create method context
32
33 Level: advanced
34
35 Note:
36 `TSGLLEAdaptRegister()` may be called multiple times to add several user-defined families.
37
38 Example Usage:
39 .vb
40 TSGLLEAdaptRegister("my_scheme", MySchemeCreate);
41 .ve
42
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
51
52 .seealso: [](ch_ts), `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegisterAll()`
53 @*/
TSGLLEAdaptRegister(const char sname[],PetscErrorCode (* function)(TSGLLEAdapt))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 }
61
62 /*@C
63 TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in `TSGLLEAdapt`
64
65 Not Collective
66
67 Level: advanced
68
69 .seealso: [](ch_ts), `TSGLLEAdapt`, `TSGLLE`, `TSGLLEAdaptRegisterDestroy()`
70 @*/
TSGLLEAdaptRegisterAll(void)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 }
81
82 /*@C
83 TSGLLEAdaptFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
84 called from `PetscFinalize()`.
85
86 Level: developer
87
88 .seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEAdapt`, `TSGLLEAdaptInitializePackage()`
89 @*/
TSGLLEAdaptFinalizePackage(void)90 PetscErrorCode TSGLLEAdaptFinalizePackage(void)
91 {
92 PetscFunctionBegin;
93 PetscCall(PetscFunctionListDestroy(&TSGLLEAdaptList));
94 TSGLLEAdaptPackageInitialized = PETSC_FALSE;
95 TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
99 /*@C
100 TSGLLEAdaptInitializePackage - This function initializes everything in the `TSGLLEAdapt` package. It is
101 called from `TSInitializePackage()`.
102
103 Level: developer
104
105 .seealso: [](ch_ts), `PetscInitialize()`, `TSGLLEAdapt`, `TSGLLEAdaptFinalizePackage()`
106 @*/
TSGLLEAdaptInitializePackage(void)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 }
117
TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)118 PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt, TSGLLEAdaptType type)
119 {
120 PetscErrorCode (*r)(TSGLLEAdapt);
121
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 }
130
TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])131 PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt, const char prefix[])
132 {
133 PetscFunctionBegin;
134 PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adapt, prefix));
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)138 PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt, PetscViewer viewer)
139 {
140 PetscBool isascii;
141
142 PetscFunctionBegin;
143 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
144 if (isascii) {
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 }
154
TSGLLEAdaptDestroy(TSGLLEAdapt * adapt)155 PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
156 {
157 PetscFunctionBegin;
158 if (!*adapt) PetscFunctionReturn(PETSC_SUCCESS);
159 PetscValidHeaderSpecific(*adapt, TSGLLEADAPT_CLASSID, 1);
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 }
168
TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt,PetscOptionItems PetscOptionsObject)169 PetscErrorCode TSGLLEAdaptSetFromOptions(TSGLLEAdapt adapt, PetscOptionItems PetscOptionsObject)
170 {
171 char type[256] = TSGLLEADAPT_BOTH;
172 PetscBool flg;
173
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 }
184
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)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;
188 PetscValidHeaderSpecific(adapt, TSGLLEADAPT_CLASSID, 1);
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 }
198
TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt * inadapt)199 PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm, TSGLLEAdapt *inadapt)
200 {
201 TSGLLEAdapt adapt;
202
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 }
209
210 /*
211 Implementations
212 */
213
TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)214 static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
215 {
216 PetscFunctionBegin;
217 PetscCall(PetscFree(adapt->data));
218 PetscFunctionReturn(PETSC_SUCCESS);
219 }
220
221 typedef struct {
222 PetscInt scheme;
223 PetscReal h;
224 } TSGLLEAdapt_None;
225
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)226 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)
227 {
228 PetscFunctionBegin;
229 *next_sc = cur;
230 *next_h = h;
231 if (*next_h > tleft) {
232 *finish = PETSC_TRUE;
233 *next_h = tleft;
234 } else *finish = PETSC_FALSE;
235 PetscFunctionReturn(PETSC_SUCCESS);
236 }
237
TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)238 PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
239 {
240 TSGLLEAdapt_None *a;
241
242 PetscFunctionBegin;
243 PetscCall(PetscNew(&a));
244 adapt->data = (void *)a;
245 adapt->ops->choose = TSGLLEAdaptChoose_None;
246 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
247 PetscFunctionReturn(PETSC_SUCCESS);
248 }
249
250 typedef struct {
251 PetscReal desired_h;
252 } TSGLLEAdapt_Size;
253
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 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)
255 {
256 TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size *)adapt->data;
257 PetscReal dec = 0.2, inc = 5.0, safe = 0.9, optimal, last_desired_h;
258
259 PetscFunctionBegin;
260 *next_sc = cur;
261 optimal = PetscPowReal((PetscReal)errors[cur], (PetscReal)-1. / (safe * orders[cur]));
262 /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
263 * one that would have been taken (without smoothing) on the last step. */
264 last_desired_h = sz->desired_h;
265 sz->desired_h = h * PetscMax(dec, PetscMin(inc, optimal)); /* Trim to [dec,inc] */
266
267 /* Normally only happens on the first step */
268 if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
269 else *next_h = sz->desired_h;
270
271 if (*next_h > tleft) {
272 *finish = PETSC_TRUE;
273 *next_h = tleft;
274 } else *finish = PETSC_FALSE;
275 PetscFunctionReturn(PETSC_SUCCESS);
276 }
277
TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)278 PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
279 {
280 TSGLLEAdapt_Size *a;
281
282 PetscFunctionBegin;
283 PetscCall(PetscNew(&a));
284 adapt->data = (void *)a;
285 adapt->ops->choose = TSGLLEAdaptChoose_Size;
286 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
287 PetscFunctionReturn(PETSC_SUCCESS);
288 }
289
290 typedef struct {
291 PetscInt count_at_order;
292 PetscReal desired_h;
293 } TSGLLEAdapt_Both;
294
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 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;
304
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 }
331
332 if (*next_h > tleft) {
333 *finish = PETSC_TRUE;
334 *next_h = tleft;
335 } else *finish = PETSC_FALSE;
336 PetscFunctionReturn(PETSC_SUCCESS);
337 }
338
TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)339 PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
340 {
341 TSGLLEAdapt_Both *a;
342
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 }
350