Actual source code: glleadapt.c
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: .seealso: TSGLLEAdaptRegisterAll()
50: @*/
51: PetscErrorCode TSGLLEAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLLEAdapt))
52: {
53: TSGLLEAdaptInitializePackage();
54: PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);
55: return 0;
56: }
58: /*@C
59: TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLLEAdapt
61: Not Collective
63: Level: advanced
65: .seealso: TSGLLEAdaptRegisterDestroy()
66: @*/
67: PetscErrorCode TSGLLEAdaptRegisterAll(void)
68: {
69: if (TSGLLEAdaptRegisterAllCalled) return 0;
70: TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
71: TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
72: TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
73: TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
74: return 0;
75: }
77: /*@C
78: TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
79: called from PetscFinalize().
81: Level: developer
83: .seealso: PetscFinalize()
84: @*/
85: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
86: {
87: PetscFunctionListDestroy(&TSGLLEAdaptList);
88: TSGLLEAdaptPackageInitialized = PETSC_FALSE;
89: TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
90: return 0;
91: }
93: /*@C
94: TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
95: called from TSInitializePackage().
97: Level: developer
99: .seealso: PetscInitialize()
100: @*/
101: PetscErrorCode TSGLLEAdaptInitializePackage(void)
102: {
103: if (TSGLLEAdaptPackageInitialized) return 0;
104: TSGLLEAdaptPackageInitialized = PETSC_TRUE;
105: PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
106: TSGLLEAdaptRegisterAll();
107: PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
108: return 0;
109: }
111: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
112: {
113: PetscErrorCode (*r)(TSGLLEAdapt);
115: PetscFunctionListFind(TSGLLEAdaptList,type,&r);
117: if (((PetscObject)adapt)->type_name) (*adapt->ops->destroy)(adapt);
118: (*r)(adapt);
119: PetscObjectChangeTypeName((PetscObject)adapt,type);
120: return 0;
121: }
123: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
124: {
125: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
126: return 0;
127: }
129: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
130: {
131: PetscBool iascii;
133: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
134: if (iascii) {
135: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
136: if (adapt->ops->view) {
137: PetscViewerASCIIPushTab(viewer);
138: (*adapt->ops->view)(adapt,viewer);
139: PetscViewerASCIIPopTab(viewer);
140: }
141: }
142: return 0;
143: }
145: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
146: {
147: if (!*adapt) return 0;
149: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return 0;}
150: if ((*adapt)->ops->destroy) (*(*adapt)->ops->destroy)(*adapt);
151: PetscHeaderDestroy(adapt);
152: return 0;
153: }
155: PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
156: {
157: char type[256] = TSGLLEADAPT_BOTH;
158: PetscBool flg;
160: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
161: * function can only be called from inside TSSetFromOptions_GLLE() */
162: PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");
163: PetscCall(PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList,
164: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg));
165: if (flg || !((PetscObject)adapt)->type_name) TSGLLEAdaptSetType(adapt,type);
166: if (adapt->ops->setfromoptions) (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);
167: PetscOptionsTail();
168: return 0;
169: }
171: 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)
172: {
180: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
181: return 0;
182: }
184: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
185: {
186: TSGLLEAdapt adapt;
188: *inadapt = NULL;
189: PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
190: *inadapt = adapt;
191: return 0;
192: }
194: /*
195: Implementations
196: */
198: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
199: {
200: PetscFree(adapt->data);
201: return 0;
202: }
204: /* -------------------------------- None ----------------------------------- */
205: typedef struct {
206: PetscInt scheme;
207: PetscReal h;
208: } TSGLLEAdapt_None;
210: 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)
211: {
212: *next_sc = cur;
213: *next_h = h;
214: if (*next_h > tleft) {
215: *finish = PETSC_TRUE;
216: *next_h = tleft;
217: } else *finish = PETSC_FALSE;
218: return 0;
219: }
221: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
222: {
223: TSGLLEAdapt_None *a;
225: PetscNewLog(adapt,&a);
226: adapt->data = (void*)a;
227: adapt->ops->choose = TSGLLEAdaptChoose_None;
228: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
229: return 0;
230: }
232: /* -------------------------------- Size ----------------------------------- */
233: typedef struct {
234: PetscReal desired_h;
235: } TSGLLEAdapt_Size;
237: 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)
238: {
239: TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
240: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
242: *next_sc = cur;
243: optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
244: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
245: * one that would have been taken (without smoothing) on the last step. */
246: last_desired_h = sz->desired_h;
247: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
249: /* Normally only happens on the first step */
250: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
251: else *next_h = sz->desired_h;
253: if (*next_h > tleft) {
254: *finish = PETSC_TRUE;
255: *next_h = tleft;
256: } else *finish = PETSC_FALSE;
257: return 0;
258: }
260: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
261: {
262: TSGLLEAdapt_Size *a;
264: PetscNewLog(adapt,&a);
265: adapt->data = (void*)a;
266: adapt->ops->choose = TSGLLEAdaptChoose_Size;
267: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
268: return 0;
269: }
271: /* -------------------------------- Both ----------------------------------- */
272: typedef struct {
273: PetscInt count_at_order;
274: PetscReal desired_h;
275: } TSGLLEAdapt_Both;
277: 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)
278: {
279: TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
280: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
281: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
282: PetscInt i;
284: for (i=0; i<n; i++) {
285: PetscReal optimal;
286: trial.id = i;
287: optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
288: trial.h = h*optimal;
289: trial.eff = trial.h/cost[i];
290: if (trial.eff > best.eff) PetscArraycpy(&best,&trial,1);
291: if (i == cur) PetscArraycpy(¤t,&trial,1);
292: }
293: /* Only switch orders if the scheme offers significant benefits over the current one.
294: When the scheme is not changing, only change step size if it offers significant benefits. */
295: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
296: PetscReal last_desired_h;
297: *next_sc = current.id;
298: last_desired_h = both->desired_h;
299: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
300: *next_h = (both->count_at_order > 0)
301: ? PetscSqrtReal(last_desired_h * both->desired_h)
302: : both->desired_h;
303: both->count_at_order++;
304: } else {
305: PetscReal rat = cost[best.id]/cost[cur];
306: *next_sc = best.id;
307: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
308: both->count_at_order = 0;
309: both->desired_h = best.h;
310: }
312: if (*next_h > tleft) {
313: *finish = PETSC_TRUE;
314: *next_h = tleft;
315: } else *finish = PETSC_FALSE;
316: return 0;
317: }
319: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
320: {
321: TSGLLEAdapt_Both *a;
323: PetscNewLog(adapt,&a);
324: adapt->data = (void*)a;
325: adapt->ops->choose = TSGLLEAdaptChoose_Both;
326: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
327: return 0;
328: }