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: {
56: TSGLLEAdaptInitializePackage();
57: PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);
58: return(0);
59: }
61: /*@C
62: TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLLEAdapt
64: Not Collective
66: Level: advanced
68: .seealso: TSGLLEAdaptRegisterDestroy()
69: @*/
70: PetscErrorCode TSGLLEAdaptRegisterAll(void)
71: {
75: if (TSGLLEAdaptRegisterAllCalled) return(0);
76: TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
77: TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
78: TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
79: TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
80: return(0);
81: }
83: /*@C
84: TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
85: called from PetscFinalize().
87: Level: developer
89: .seealso: PetscFinalize()
90: @*/
91: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
92: {
96: PetscFunctionListDestroy(&TSGLLEAdaptList);
97: TSGLLEAdaptPackageInitialized = PETSC_FALSE;
98: TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
99: return(0);
100: }
102: /*@C
103: TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
104: called from TSInitializePackage().
106: Level: developer
108: .seealso: PetscInitialize()
109: @*/
110: PetscErrorCode TSGLLEAdaptInitializePackage(void)
111: {
115: if (TSGLLEAdaptPackageInitialized) return(0);
116: TSGLLEAdaptPackageInitialized = PETSC_TRUE;
117: PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
118: TSGLLEAdaptRegisterAll();
119: PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
120: return(0);
121: }
123: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
124: {
125: PetscErrorCode ierr,(*r)(TSGLLEAdapt);
128: PetscFunctionListFind(TSGLLEAdaptList,type,&r);
129: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type);
130: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
131: (*r)(adapt);
132: PetscObjectChangeTypeName((PetscObject)adapt,type);
133: return(0);
134: }
136: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
137: {
141: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
142: return(0);
143: }
145: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
146: {
148: PetscBool iascii;
151: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
152: if (iascii) {
153: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
154: if (adapt->ops->view) {
155: PetscViewerASCIIPushTab(viewer);
156: (*adapt->ops->view)(adapt,viewer);
157: PetscViewerASCIIPopTab(viewer);
158: }
159: }
160: return(0);
161: }
163: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
164: {
168: if (!*adapt) return(0);
170: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; return(0);}
171: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
172: PetscHeaderDestroy(adapt);
173: return(0);
174: }
176: PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
177: {
179: char type[256] = TSGLLEADAPT_BOTH;
180: PetscBool flg;
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: PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");
186: 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);
188: if (flg || !((PetscObject)adapt)->type_name) {
189: TSGLLEAdaptSetType(adapt,type);
190: }
191: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
192: PetscOptionsTail();
193: return(0);
194: }
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: {
208: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
209: return(0);
210: }
212: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
213: {
215: TSGLLEAdapt adapt;
218: *inadapt = NULL;
219: PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
220: *inadapt = adapt;
221: return(0);
222: }
224: /*
225: Implementations
226: */
228: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
229: {
233: PetscFree(adapt->data);
234: return(0);
235: }
237: /* -------------------------------- None ----------------------------------- */
238: typedef struct {
239: PetscInt scheme;
240: PetscReal h;
241: } TSGLLEAdapt_None;
243: 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)
244: {
247: *next_sc = cur;
248: *next_h = h;
249: if (*next_h > tleft) {
250: *finish = PETSC_TRUE;
251: *next_h = tleft;
252: } else *finish = PETSC_FALSE;
253: return(0);
254: }
256: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
257: {
259: TSGLLEAdapt_None *a;
262: PetscNewLog(adapt,&a);
263: adapt->data = (void*)a;
264: adapt->ops->choose = TSGLLEAdaptChoose_None;
265: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
266: return(0);
267: }
269: /* -------------------------------- Size ----------------------------------- */
270: typedef struct {
271: PetscReal desired_h;
272: } TSGLLEAdapt_Size;
274: 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)
275: {
276: TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
277: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
280: *next_sc = cur;
281: optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
282: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
283: * one that would have been taken (without smoothing) on the last step. */
284: last_desired_h = sz->desired_h;
285: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
287: /* Normally only happens on the first step */
288: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
289: else *next_h = sz->desired_h;
291: if (*next_h > tleft) {
292: *finish = PETSC_TRUE;
293: *next_h = tleft;
294: } else *finish = PETSC_FALSE;
295: return(0);
296: }
298: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
299: {
301: TSGLLEAdapt_Size *a;
304: PetscNewLog(adapt,&a);
305: adapt->data = (void*)a;
306: adapt->ops->choose = TSGLLEAdaptChoose_Size;
307: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
308: return(0);
309: }
311: /* -------------------------------- Both ----------------------------------- */
312: typedef struct {
313: PetscInt count_at_order;
314: PetscReal desired_h;
315: } TSGLLEAdapt_Both;
317: 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)
318: {
319: TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
320: PetscErrorCode ierr;
321: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
322: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
323: PetscInt i;
326: for (i=0; i<n; i++) {
327: PetscReal optimal;
328: trial.id = i;
329: optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
330: trial.h = h*optimal;
331: trial.eff = trial.h/cost[i];
332: if (trial.eff > best.eff) {PetscArraycpy(&best,&trial,1);}
333: if (i == cur) {PetscArraycpy(¤t,&trial,1);}
334: }
335: /* Only switch orders if the scheme offers significant benefits over the current one.
336: When the scheme is not changing, only change step size if it offers significant benefits. */
337: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
338: PetscReal last_desired_h;
339: *next_sc = current.id;
340: last_desired_h = both->desired_h;
341: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
342: *next_h = (both->count_at_order > 0)
343: ? PetscSqrtReal(last_desired_h * both->desired_h)
344: : both->desired_h;
345: both->count_at_order++;
346: } else {
347: PetscReal rat = cost[best.id]/cost[cur];
348: *next_sc = best.id;
349: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
350: both->count_at_order = 0;
351: both->desired_h = best.h;
352: }
354: if (*next_h > tleft) {
355: *finish = PETSC_TRUE;
356: *next_h = tleft;
357: } else *finish = PETSC_FALSE;
358: return(0);
359: }
361: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
362: {
364: TSGLLEAdapt_Both *a;
367: PetscNewLog(adapt,&a);
368: adapt->data = (void*)a;
369: adapt->ops->choose = TSGLLEAdaptChoose_Both;
370: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
371: return(0);
372: }