Actual source code: gladapt.c
petsc-3.5.4 2015-05-23
2: #include <../src/ts/impls/implicit/gl/gl.h> /*I "petscts.h" I*/
4: static PetscFunctionList TSGLAdaptList;
5: static PetscBool TSGLAdaptPackageInitialized;
6: static PetscBool TSGLAdaptRegisterAllCalled;
7: static PetscClassId TSGLADAPT_CLASSID;
9: struct _TSGLAdaptOps {
10: PetscErrorCode (*choose)(TSGLAdapt,PetscInt,const PetscInt[],const PetscReal[],const PetscReal[],PetscInt,PetscReal,PetscReal,PetscInt*,PetscReal*,PetscBool*);
11: PetscErrorCode (*destroy)(TSGLAdapt);
12: PetscErrorCode (*view)(TSGLAdapt,PetscViewer);
13: PetscErrorCode (*setfromoptions)(TSGLAdapt);
14: };
16: struct _p_TSGLAdapt {
17: PETSCHEADER(struct _TSGLAdaptOps);
18: void *data;
19: };
21: PETSC_EXTERN PetscErrorCode TSGLAdaptCreate_None(TSGLAdapt);
22: PETSC_EXTERN PetscErrorCode TSGLAdaptCreate_Size(TSGLAdapt);
23: PETSC_EXTERN PetscErrorCode TSGLAdaptCreate_Both(TSGLAdapt);
27: /*@C
28: TSGLAdaptRegister - adds a TSGLAdapt implementation
30: Not Collective
32: Input Parameters:
33: + name_scheme - name of user-defined adaptivity scheme
34: - routine_create - routine to create method context
36: Notes:
37: TSGLAdaptRegister() may be called multiple times to add several user-defined families.
39: Sample usage:
40: .vb
41: TSGLAdaptRegister("my_scheme",MySchemeCreate);
42: .ve
44: Then, your scheme can be chosen with the procedural interface via
45: $ TSGLAdaptSetType(ts,"my_scheme")
46: or at runtime via the option
47: $ -ts_adapt_type my_scheme
49: Level: advanced
51: .keywords: TSGLAdapt, register
53: .seealso: TSGLAdaptRegisterAll()
54: @*/
55: PetscErrorCode TSGLAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLAdapt))
56: {
60: PetscFunctionListAdd(&TSGLAdaptList,sname,function);
61: return(0);
62: }
66: /*@C
67: TSGLAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLAdapt
69: Not Collective
71: Level: advanced
73: .keywords: TSGLAdapt, register, all
75: .seealso: TSGLAdaptRegisterDestroy()
76: @*/
77: PetscErrorCode TSGLAdaptRegisterAll(void)
78: {
82: TSGLAdaptRegister(TSGLADAPT_NONE,TSGLAdaptCreate_None);
83: TSGLAdaptRegister(TSGLADAPT_SIZE,TSGLAdaptCreate_Size);
84: TSGLAdaptRegister(TSGLADAPT_BOTH,TSGLAdaptCreate_Both);
85: return(0);
86: }
90: /*@C
91: TSGLFinalizePackage - This function destroys everything in the TSGL package. It is
92: called from PetscFinalize().
94: Level: developer
96: .keywords: Petsc, destroy, package
97: .seealso: PetscFinalize()
98: @*/
99: PetscErrorCode TSGLAdaptFinalizePackage(void)
100: {
104: PetscFunctionListDestroy(&TSGLAdaptList);
105: TSGLAdaptPackageInitialized = PETSC_FALSE;
106: TSGLAdaptRegisterAllCalled = PETSC_FALSE;
107: return(0);
108: }
112: /*@C
113: TSGLAdaptInitializePackage - This function initializes everything in the TSGLAdapt package. It is
114: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
115: TSCreate_GL() when using static libraries.
117: Level: developer
119: .keywords: TSGLAdapt, initialize, package
120: .seealso: PetscInitialize()
121: @*/
122: PetscErrorCode TSGLAdaptInitializePackage(void)
123: {
127: if (TSGLAdaptPackageInitialized) return(0);
128: TSGLAdaptPackageInitialized = PETSC_TRUE;
129: PetscClassIdRegister("TSGLAdapt",&TSGLADAPT_CLASSID);
130: TSGLAdaptRegisterAll();
131: PetscRegisterFinalize(TSGLAdaptFinalizePackage);
132: return(0);
133: }
137: PetscErrorCode TSGLAdaptSetType(TSGLAdapt adapt,TSGLAdaptType type)
138: {
139: PetscErrorCode ierr,(*r)(TSGLAdapt);
142: PetscFunctionListFind(TSGLAdaptList,type,&r);
143: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLAdapt type \"%s\" given",type);
144: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
145: (*r)(adapt);
146: PetscObjectChangeTypeName((PetscObject)adapt,type);
147: return(0);
148: }
152: PetscErrorCode TSGLAdaptSetOptionsPrefix(TSGLAdapt adapt,const char prefix[])
153: {
157: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
158: return(0);
159: }
163: PetscErrorCode TSGLAdaptView(TSGLAdapt adapt,PetscViewer viewer)
164: {
166: PetscBool iascii;
169: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
170: if (iascii) {
171: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
172: if (adapt->ops->view) {
173: PetscViewerASCIIPushTab(viewer);
174: (*adapt->ops->view)(adapt,viewer);
175: PetscViewerASCIIPopTab(viewer);
176: }
177: }
178: return(0);
179: }
183: PetscErrorCode TSGLAdaptDestroy(TSGLAdapt *adapt)
184: {
188: if (!*adapt) return(0);
190: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
191: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
192: PetscHeaderDestroy(adapt);
193: return(0);
194: }
198: PetscErrorCode TSGLAdaptSetFromOptions(TSGLAdapt adapt)
199: {
201: char type[256] = TSGLADAPT_BOTH;
202: PetscBool flg;
205: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGL, but currently this
206: * function can only be called from inside TSSetFromOptions_GL() */
207: PetscOptionsHead("TSGL Adaptivity options");
208: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLAdaptSetType",TSGLAdaptList,
209: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
210: if (flg || !((PetscObject)adapt)->type_name) {
211: TSGLAdaptSetType(adapt,type);
212: }
213: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(adapt);}
214: PetscOptionsTail();
215: return(0);
216: }
220: PetscErrorCode TSGLAdaptChoose(TSGLAdapt 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)
221: {
232: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
233: return(0);
234: }
238: PetscErrorCode TSGLAdaptCreate(MPI_Comm comm,TSGLAdapt *inadapt)
239: {
241: TSGLAdapt adapt;
244: *inadapt = 0;
245: PetscHeaderCreate(adapt,_p_TSGLAdapt,struct _TSGLAdaptOps,TSGLADAPT_CLASSID,"TSGLAdapt","General Linear adaptivity","TS",comm,TSGLAdaptDestroy,TSGLAdaptView);
246: *inadapt = adapt;
247: return(0);
248: }
251: /*
252: * Implementations
253: */
257: static PetscErrorCode TSGLAdaptDestroy_JustFree(TSGLAdapt adapt)
258: {
262: PetscFree(adapt->data);
263: return(0);
264: }
266: /* -------------------------------- None ----------------------------------- */
267: typedef struct {
268: PetscInt scheme;
269: PetscReal h;
270: } TSGLAdapt_None;
274: static PetscErrorCode TSGLAdaptChoose_None(TSGLAdapt 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: {
278: *next_sc = cur;
279: *next_h = h;
280: if (*next_h > tleft) {
281: *finish = PETSC_TRUE;
282: *next_h = tleft;
283: } else *finish = PETSC_FALSE;
284: return(0);
285: }
289: PetscErrorCode TSGLAdaptCreate_None(TSGLAdapt adapt)
290: {
292: TSGLAdapt_None *a;
295: PetscNewLog(adapt,&a);
296: adapt->data = (void*)a;
297: adapt->ops->choose = TSGLAdaptChoose_None;
298: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
299: return(0);
300: }
302: /* -------------------------------- Size ----------------------------------- */
303: typedef struct {
304: PetscReal desired_h;
305: } TSGLAdapt_Size;
310: static PetscErrorCode TSGLAdaptChoose_Size(TSGLAdapt 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)
311: {
312: TSGLAdapt_Size *sz = (TSGLAdapt_Size*)adapt->data;
313: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
316: *next_sc = cur;
317: optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
318: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
319: * one that would have been taken (without smoothing) on the last step. */
320: last_desired_h = sz->desired_h;
321: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
323: /* Normally only happens on the first step */
324: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
325: else *next_h = sz->desired_h;
327: if (*next_h > tleft) {
328: *finish = PETSC_TRUE;
329: *next_h = tleft;
330: } else *finish = PETSC_FALSE;
331: return(0);
332: }
336: PetscErrorCode TSGLAdaptCreate_Size(TSGLAdapt adapt)
337: {
339: TSGLAdapt_Size *a;
342: PetscNewLog(adapt,&a);
343: adapt->data = (void*)a;
344: adapt->ops->choose = TSGLAdaptChoose_Size;
345: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
346: return(0);
347: }
349: /* -------------------------------- Both ----------------------------------- */
350: typedef struct {
351: PetscInt count_at_order;
352: PetscReal desired_h;
353: } TSGLAdapt_Both;
358: static PetscErrorCode TSGLAdaptChoose_Both(TSGLAdapt 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)
359: {
360: TSGLAdapt_Both *both = (TSGLAdapt_Both*)adapt->data;
362: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
363: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
364: PetscInt i;
367: for (i=0; i<n; i++) {
368: PetscReal optimal;
369: trial.id = i;
370: optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
371: trial.h = h*optimal;
372: trial.eff = trial.h/cost[i];
373: if (trial.eff > best.eff) {PetscMemcpy(&best,&trial,sizeof(trial));}
374: if (i == cur) {PetscMemcpy(¤t,&trial,sizeof(trial));}
375: }
376: /* Only switch orders if the scheme offers significant benefits over the current one.
377: When the scheme is not changing, only change step size if it offers significant benefits. */
378: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
379: PetscReal last_desired_h;
380: *next_sc = current.id;
381: last_desired_h = both->desired_h;
382: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
383: *next_h = (both->count_at_order > 0)
384: ? PetscSqrtReal(last_desired_h * both->desired_h)
385: : both->desired_h;
386: both->count_at_order++;
387: } else {
388: PetscReal rat = cost[best.id]/cost[cur];
389: *next_sc = best.id;
390: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
391: both->count_at_order = 0;
392: both->desired_h = best.h;
393: }
395: if (*next_h > tleft) {
396: *finish = PETSC_TRUE;
397: *next_h = tleft;
398: } else *finish = PETSC_FALSE;
399: return(0);
400: }
404: PetscErrorCode TSGLAdaptCreate_Both(TSGLAdapt adapt)
405: {
407: TSGLAdapt_Both *a;
410: PetscNewLog(adapt,&a);
411: adapt->data = (void*)a;
412: adapt->ops->choose = TSGLAdaptChoose_Both;
413: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
414: return(0);
415: }