1
2 #include <../src/ts/impls/implicit/glle/glle.h> /*I "petscts.h" I*/
3
4 static PetscFunctionList TSGLLEAdaptList;
5 static PetscBool TSGLLEAdaptPackageInitialized;
6 static PetscBool TSGLLEAdaptRegisterAllCalled;
7 static PetscClassId TSGLLEADAPT_CLASSID;
8
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)(PetscOptionItems*,TSGLLEAdapt);
14 };
15
16 struct _p_TSGLLEAdapt {
17 PETSCHEADER(struct _TSGLLEAdaptOps);
18 void *data;
19 };
20
21 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt);
22 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt);
23 PETSC_EXTERN PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt);
24
25 /*@C
26 TSGLLEAdaptRegister - adds a TSGLLEAdapt implementation
27
28 Not Collective
29
30 Input Parameters:
31 + name_scheme - name of user-defined adaptivity scheme
32 - routine_create - routine to create method context
33
34 Notes:
35 TSGLLEAdaptRegister() may be called multiple times to add several user-defined families.
36
37 Sample usage:
38 .vb
39 TSGLLEAdaptRegister("my_scheme",MySchemeCreate);
40 .ve
41
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
46
47 Level: advanced
48
49 .seealso: TSGLLEAdaptRegisterAll()
50 @*/
TSGLLEAdaptRegister(const char sname[],PetscErrorCode (* function)(TSGLLEAdapt))51 PetscErrorCode TSGLLEAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLLEAdapt))
52 {
53 PetscErrorCode ierr;
54
55 PetscFunctionBegin;
56 ierr = TSGLLEAdaptInitializePackage();CHKERRQ(ierr);
57 ierr = PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);CHKERRQ(ierr);
58 PetscFunctionReturn(0);
59 }
60
61 /*@C
62 TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLLEAdapt
63
64 Not Collective
65
66 Level: advanced
67
68 .seealso: TSGLLEAdaptRegisterDestroy()
69 @*/
TSGLLEAdaptRegisterAll(void)70 PetscErrorCode TSGLLEAdaptRegisterAll(void)
71 {
72 PetscErrorCode ierr;
73
74 PetscFunctionBegin;
75 if (TSGLLEAdaptRegisterAllCalled) PetscFunctionReturn(0);
76 TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
77 ierr = TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);CHKERRQ(ierr);
78 ierr = TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);CHKERRQ(ierr);
79 ierr = TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);CHKERRQ(ierr);
80 PetscFunctionReturn(0);
81 }
82
83 /*@C
84 TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
85 called from PetscFinalize().
86
87 Level: developer
88
89 .seealso: PetscFinalize()
90 @*/
TSGLLEAdaptFinalizePackage(void)91 PetscErrorCode TSGLLEAdaptFinalizePackage(void)
92 {
93 PetscErrorCode ierr;
94
95 PetscFunctionBegin;
96 ierr = PetscFunctionListDestroy(&TSGLLEAdaptList);CHKERRQ(ierr);
97 TSGLLEAdaptPackageInitialized = PETSC_FALSE;
98 TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
99 PetscFunctionReturn(0);
100 }
101
102 /*@C
103 TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
104 called from TSInitializePackage().
105
106 Level: developer
107
108 .seealso: PetscInitialize()
109 @*/
TSGLLEAdaptInitializePackage(void)110 PetscErrorCode TSGLLEAdaptInitializePackage(void)
111 {
112 PetscErrorCode ierr;
113
114 PetscFunctionBegin;
115 if (TSGLLEAdaptPackageInitialized) PetscFunctionReturn(0);
116 TSGLLEAdaptPackageInitialized = PETSC_TRUE;
117 ierr = PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);CHKERRQ(ierr);
118 ierr = TSGLLEAdaptRegisterAll();CHKERRQ(ierr);
119 ierr = PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);CHKERRQ(ierr);
120 PetscFunctionReturn(0);
121 }
122
TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)123 PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
124 {
125 PetscErrorCode ierr,(*r)(TSGLLEAdapt);
126
127 PetscFunctionBegin;
128 ierr = PetscFunctionListFind(TSGLLEAdaptList,type,&r);CHKERRQ(ierr);
129 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type);
130 if (((PetscObject)adapt)->type_name) {ierr = (*adapt->ops->destroy)(adapt);CHKERRQ(ierr);}
131 ierr = (*r)(adapt);CHKERRQ(ierr);
132 ierr = PetscObjectChangeTypeName((PetscObject)adapt,type);CHKERRQ(ierr);
133 PetscFunctionReturn(0);
134 }
135
TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])136 PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
137 {
138 PetscErrorCode ierr;
139
140 PetscFunctionBegin;
141 ierr = PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);CHKERRQ(ierr);
142 PetscFunctionReturn(0);
143 }
144
TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)145 PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
146 {
147 PetscErrorCode ierr;
148 PetscBool iascii;
149
150 PetscFunctionBegin;
151 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
152 if (iascii) {
153 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);CHKERRQ(ierr);
154 if (adapt->ops->view) {
155 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
156 ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr);
157 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
158 }
159 }
160 PetscFunctionReturn(0);
161 }
162
TSGLLEAdaptDestroy(TSGLLEAdapt * adapt)163 PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
164 {
165 PetscErrorCode ierr;
166
167 PetscFunctionBegin;
168 if (!*adapt) PetscFunctionReturn(0);
169 PetscValidHeaderSpecific(*adapt,TSGLLEADAPT_CLASSID,1);
170 if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; PetscFunctionReturn(0);}
171 if ((*adapt)->ops->destroy) {ierr = (*(*adapt)->ops->destroy)(*adapt);CHKERRQ(ierr);}
172 ierr = PetscHeaderDestroy(adapt);CHKERRQ(ierr);
173 PetscFunctionReturn(0);
174 }
175
TSGLLEAdaptSetFromOptions(PetscOptionItems * PetscOptionsObject,TSGLLEAdapt adapt)176 PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
177 {
178 PetscErrorCode ierr;
179 char type[256] = TSGLLEADAPT_BOTH;
180 PetscBool flg;
181
182 PetscFunctionBegin;
183 /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
184 * function can only be called from inside TSSetFromOptions_GLLE() */
185 ierr = PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");CHKERRQ(ierr);
186 ierr = PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList,
187 ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);CHKERRQ(ierr);
188 if (flg || !((PetscObject)adapt)->type_name) {
189 ierr = TSGLLEAdaptSetType(adapt,type);CHKERRQ(ierr);
190 }
191 if (adapt->ops->setfromoptions) {ierr = (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);CHKERRQ(ierr);}
192 ierr = PetscOptionsTail();CHKERRQ(ierr);
193 PetscFunctionReturn(0);
194 }
195
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)196 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)
197 {
198 PetscErrorCode ierr;
199
200 PetscFunctionBegin;
201 PetscValidHeaderSpecific(adapt,TSGLLEADAPT_CLASSID,1);
202 PetscValidIntPointer(orders,3);
203 PetscValidPointer(errors,4);
204 PetscValidPointer(cost,5);
205 PetscValidIntPointer(next_sc,9);
206 PetscValidPointer(next_h,10);
207 PetscValidIntPointer(finish,11);
208 ierr = (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);CHKERRQ(ierr);
209 PetscFunctionReturn(0);
210 }
211
TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt * inadapt)212 PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
213 {
214 PetscErrorCode ierr;
215 TSGLLEAdapt adapt;
216
217 PetscFunctionBegin;
218 *inadapt = NULL;
219 ierr = PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);CHKERRQ(ierr);
220 *inadapt = adapt;
221 PetscFunctionReturn(0);
222 }
223
224
225 /*
226 * Implementations
227 */
228
TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)229 static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
230 {
231 PetscErrorCode ierr;
232
233 PetscFunctionBegin;
234 ierr = PetscFree(adapt->data);CHKERRQ(ierr);
235 PetscFunctionReturn(0);
236 }
237
238 /* -------------------------------- None ----------------------------------- */
239 typedef struct {
240 PetscInt scheme;
241 PetscReal h;
242 } TSGLLEAdapt_None;
243
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)244 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)
245 {
246
247 PetscFunctionBegin;
248 *next_sc = cur;
249 *next_h = h;
250 if (*next_h > tleft) {
251 *finish = PETSC_TRUE;
252 *next_h = tleft;
253 } else *finish = PETSC_FALSE;
254 PetscFunctionReturn(0);
255 }
256
TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)257 PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
258 {
259 PetscErrorCode ierr;
260 TSGLLEAdapt_None *a;
261
262 PetscFunctionBegin;
263 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr);
264 adapt->data = (void*)a;
265 adapt->ops->choose = TSGLLEAdaptChoose_None;
266 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
267 PetscFunctionReturn(0);
268 }
269
270 /* -------------------------------- Size ----------------------------------- */
271 typedef struct {
272 PetscReal desired_h;
273 } TSGLLEAdapt_Size;
274
275
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)276 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)
277 {
278 TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
279 PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
280
281 PetscFunctionBegin;
282 *next_sc = cur;
283 optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
284 /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
285 * one that would have been taken (without smoothing) on the last step. */
286 last_desired_h = sz->desired_h;
287 sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
288
289 /* Normally only happens on the first step */
290 if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
291 else *next_h = sz->desired_h;
292
293 if (*next_h > tleft) {
294 *finish = PETSC_TRUE;
295 *next_h = tleft;
296 } else *finish = PETSC_FALSE;
297 PetscFunctionReturn(0);
298 }
299
TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)300 PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
301 {
302 PetscErrorCode ierr;
303 TSGLLEAdapt_Size *a;
304
305 PetscFunctionBegin;
306 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr);
307 adapt->data = (void*)a;
308 adapt->ops->choose = TSGLLEAdaptChoose_Size;
309 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
310 PetscFunctionReturn(0);
311 }
312
313 /* -------------------------------- Both ----------------------------------- */
314 typedef struct {
315 PetscInt count_at_order;
316 PetscReal desired_h;
317 } TSGLLEAdapt_Both;
318
319
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)320 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)
321 {
322 TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
323 PetscErrorCode ierr;
324 PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
325 struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
326 PetscInt i;
327
328 PetscFunctionBegin;
329 for (i=0; i<n; i++) {
330 PetscReal optimal;
331 trial.id = i;
332 optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
333 trial.h = h*optimal;
334 trial.eff = trial.h/cost[i];
335 if (trial.eff > best.eff) {ierr = PetscArraycpy(&best,&trial,1);CHKERRQ(ierr);}
336 if (i == cur) {ierr = PetscArraycpy(¤t,&trial,1);CHKERRQ(ierr);}
337 }
338 /* Only switch orders if the scheme offers significant benefits over the current one.
339 When the scheme is not changing, only change step size if it offers significant benefits. */
340 if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
341 PetscReal last_desired_h;
342 *next_sc = current.id;
343 last_desired_h = both->desired_h;
344 both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
345 *next_h = (both->count_at_order > 0)
346 ? PetscSqrtReal(last_desired_h * both->desired_h)
347 : both->desired_h;
348 both->count_at_order++;
349 } else {
350 PetscReal rat = cost[best.id]/cost[cur];
351 *next_sc = best.id;
352 *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
353 both->count_at_order = 0;
354 both->desired_h = best.h;
355 }
356
357 if (*next_h > tleft) {
358 *finish = PETSC_TRUE;
359 *next_h = tleft;
360 } else *finish = PETSC_FALSE;
361 PetscFunctionReturn(0);
362 }
363
TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)364 PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
365 {
366 PetscErrorCode ierr;
367 TSGLLEAdapt_Both *a;
368
369 PetscFunctionBegin;
370 ierr = PetscNewLog(adapt,&a);CHKERRQ(ierr);
371 adapt->data = (void*)a;
372 adapt->ops->choose = TSGLLEAdaptChoose_Both;
373 adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
374 PetscFunctionReturn(0);
375 }
376