Actual source code: glleadapt.c
petsc-3.11.4 2019-09-28
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)(PetscOptionItems*,TSGLLEAdapt);
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: + name_scheme - name of user-defined adaptivity scheme
32: - routine_create - routine to create method context
34: Notes:
35: TSGLLEAdaptRegister() may be called multiple times to add several user-defined families.
37: Sample usage:
38: .vb
39: TSGLLEAdaptRegister("my_scheme",MySchemeCreate);
40: .ve
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
47: Level: advanced
49: .keywords: TSGLLEAdapt, register
51: .seealso: TSGLLEAdaptRegisterAll()
52: @*/
53: PetscErrorCode TSGLLEAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLLEAdapt))
54: {
58: TSGLLEAdaptInitializePackage();
59: PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);
60: return(0);
61: }
63: /*@C
64: TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLLEAdapt
66: Not Collective
68: Level: advanced
70: .keywords: TSGLLEAdapt, register, all
72: .seealso: TSGLLEAdaptRegisterDestroy()
73: @*/
74: PetscErrorCode TSGLLEAdaptRegisterAll(void)
75: {
79: if (TSGLLEAdaptRegisterAllCalled) return(0);
80: TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
81: TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
82: TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
83: TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
84: return(0);
85: }
87: /*@C
88: TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
89: called from PetscFinalize().
91: Level: developer
93: .keywords: Petsc, destroy, package
94: .seealso: PetscFinalize()
95: @*/
96: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
97: {
101: PetscFunctionListDestroy(&TSGLLEAdaptList);
102: TSGLLEAdaptPackageInitialized = PETSC_FALSE;
103: TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
104: return(0);
105: }
107: /*@C
108: TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
109: called from TSInitializePackage().
111: Level: developer
113: .keywords: TSGLLEAdapt, initialize, package
114: .seealso: PetscInitialize()
115: @*/
116: PetscErrorCode TSGLLEAdaptInitializePackage(void)
117: {
121: if (TSGLLEAdaptPackageInitialized) return(0);
122: TSGLLEAdaptPackageInitialized = PETSC_TRUE;
123: PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
124: TSGLLEAdaptRegisterAll();
125: PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
126: return(0);
127: }
129: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
130: {
131: PetscErrorCode ierr,(*r)(TSGLLEAdapt);
134: PetscFunctionListFind(TSGLLEAdaptList,type,&r);
135: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type);
136: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
137: (*r)(adapt);
138: PetscObjectChangeTypeName((PetscObject)adapt,type);
139: return(0);
140: }
142: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
143: {
147: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
148: return(0);
149: }
151: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
152: {
154: PetscBool iascii;
157: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
158: if (iascii) {
159: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
160: if (adapt->ops->view) {
161: PetscViewerASCIIPushTab(viewer);
162: (*adapt->ops->view)(adapt,viewer);
163: PetscViewerASCIIPopTab(viewer);
164: }
165: }
166: return(0);
167: }
169: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
170: {
174: if (!*adapt) return(0);
176: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
177: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
178: PetscHeaderDestroy(adapt);
179: return(0);
180: }
182: PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
183: {
185: char type[256] = TSGLLEADAPT_BOTH;
186: PetscBool flg;
189: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
190: * function can only be called from inside TSSetFromOptions_GLLE() */
191: PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");
192: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList,
193: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
194: if (flg || !((PetscObject)adapt)->type_name) {
195: TSGLLEAdaptSetType(adapt,type);
196: }
197: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
198: PetscOptionsTail();
199: return(0);
200: }
202: 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)
203: {
214: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
215: return(0);
216: }
218: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
219: {
221: TSGLLEAdapt adapt;
224: *inadapt = NULL;
225: PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
226: *inadapt = adapt;
227: return(0);
228: }
231: /*
232: * Implementations
233: */
235: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
236: {
240: PetscFree(adapt->data);
241: return(0);
242: }
244: /* -------------------------------- None ----------------------------------- */
245: typedef struct {
246: PetscInt scheme;
247: PetscReal h;
248: } TSGLLEAdapt_None;
250: 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)
251: {
254: *next_sc = cur;
255: *next_h = h;
256: if (*next_h > tleft) {
257: *finish = PETSC_TRUE;
258: *next_h = tleft;
259: } else *finish = PETSC_FALSE;
260: return(0);
261: }
263: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
264: {
266: TSGLLEAdapt_None *a;
269: PetscNewLog(adapt,&a);
270: adapt->data = (void*)a;
271: adapt->ops->choose = TSGLLEAdaptChoose_None;
272: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
273: return(0);
274: }
276: /* -------------------------------- Size ----------------------------------- */
277: typedef struct {
278: PetscReal desired_h;
279: } TSGLLEAdapt_Size;
282: 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)
283: {
284: TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
285: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
288: *next_sc = cur;
289: optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
290: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
291: * one that would have been taken (without smoothing) on the last step. */
292: last_desired_h = sz->desired_h;
293: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
295: /* Normally only happens on the first step */
296: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
297: else *next_h = sz->desired_h;
299: if (*next_h > tleft) {
300: *finish = PETSC_TRUE;
301: *next_h = tleft;
302: } else *finish = PETSC_FALSE;
303: return(0);
304: }
306: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
307: {
309: TSGLLEAdapt_Size *a;
312: PetscNewLog(adapt,&a);
313: adapt->data = (void*)a;
314: adapt->ops->choose = TSGLLEAdaptChoose_Size;
315: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
316: return(0);
317: }
319: /* -------------------------------- Both ----------------------------------- */
320: typedef struct {
321: PetscInt count_at_order;
322: PetscReal desired_h;
323: } TSGLLEAdapt_Both;
326: 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)
327: {
328: TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
330: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
331: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
332: PetscInt i;
335: for (i=0; i<n; i++) {
336: PetscReal optimal;
337: trial.id = i;
338: optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
339: trial.h = h*optimal;
340: trial.eff = trial.h/cost[i];
341: if (trial.eff > best.eff) {PetscMemcpy(&best,&trial,sizeof(trial));}
342: if (i == cur) {PetscMemcpy(¤t,&trial,sizeof(trial));}
343: }
344: /* Only switch orders if the scheme offers significant benefits over the current one.
345: When the scheme is not changing, only change step size if it offers significant benefits. */
346: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
347: PetscReal last_desired_h;
348: *next_sc = current.id;
349: last_desired_h = both->desired_h;
350: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
351: *next_h = (both->count_at_order > 0)
352: ? PetscSqrtReal(last_desired_h * both->desired_h)
353: : both->desired_h;
354: both->count_at_order++;
355: } else {
356: PetscReal rat = cost[best.id]/cost[cur];
357: *next_sc = best.id;
358: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
359: both->count_at_order = 0;
360: both->desired_h = best.h;
361: }
363: if (*next_h > tleft) {
364: *finish = PETSC_TRUE;
365: *next_h = tleft;
366: } else *finish = PETSC_FALSE;
367: return(0);
368: }
370: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
371: {
373: TSGLLEAdapt_Both *a;
376: PetscNewLog(adapt,&a);
377: adapt->data = (void*)a;
378: adapt->ops->choose = TSGLLEAdaptChoose_Both;
379: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
380: return(0);
381: }