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: }
225: /*
226: * Implementations
227: */
229: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
230: {
234: PetscFree(adapt->data);
235: return(0);
236: }
238: /* -------------------------------- None ----------------------------------- */
239: typedef struct {
240: PetscInt scheme;
241: PetscReal h;
242: } TSGLLEAdapt_None;
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: {
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: return(0);
255: }
257: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
258: {
260: TSGLLEAdapt_None *a;
263: PetscNewLog(adapt,&a);
264: adapt->data = (void*)a;
265: adapt->ops->choose = TSGLLEAdaptChoose_None;
266: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
267: return(0);
268: }
270: /* -------------------------------- Size ----------------------------------- */
271: typedef struct {
272: PetscReal desired_h;
273: } TSGLLEAdapt_Size;
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;
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] */
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;
293: if (*next_h > tleft) {
294: *finish = PETSC_TRUE;
295: *next_h = tleft;
296: } else *finish = PETSC_FALSE;
297: return(0);
298: }
300: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
301: {
303: TSGLLEAdapt_Size *a;
306: PetscNewLog(adapt,&a);
307: adapt->data = (void*)a;
308: adapt->ops->choose = TSGLLEAdaptChoose_Size;
309: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
310: return(0);
311: }
313: /* -------------------------------- Both ----------------------------------- */
314: typedef struct {
315: PetscInt count_at_order;
316: PetscReal desired_h;
317: } TSGLLEAdapt_Both;
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;
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) {PetscArraycpy(&best,&trial,1);}
336: if (i == cur) {PetscArraycpy(¤t,&trial,1);}
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: }
357: if (*next_h > tleft) {
358: *finish = PETSC_TRUE;
359: *next_h = tleft;
360: } else *finish = PETSC_FALSE;
361: return(0);
362: }
364: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
365: {
367: TSGLLEAdapt_Both *a;
370: PetscNewLog(adapt,&a);
371: adapt->data = (void*)a;
372: adapt->ops->choose = TSGLLEAdaptChoose_Both;
373: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
374: return(0);
375: }