Actual source code: gladapt.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ts/impls/implicit/gl/gl.h> /*I "petscts.h" I*/
4: static PetscFList 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: EXTERN_C_BEGIN
22: PetscErrorCode TSGLAdaptCreate_None(TSGLAdapt);
23: PetscErrorCode TSGLAdaptCreate_Size(TSGLAdapt);
24: PetscErrorCode TSGLAdaptCreate_Both(TSGLAdapt);
25: EXTERN_C_END
29: /*@C
30: TSGLAdaptRegister - see TSGLAdaptRegisterDynamic()
32: Level: advanced
33: @*/
34: PetscErrorCode TSGLAdaptRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(TSGLAdapt))
35: {
37: char fullname[PETSC_MAX_PATH_LEN];
40: PetscFListConcat(path,name,fullname);
41: PetscFListAdd(&TSGLAdaptList,sname,fullname,(void(*)(void))function);
42: return(0);
43: }
47: /*@C
48: TSGLAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLAdapt
50: Not Collective
52: Level: advanced
54: .keywords: TSGLAdapt, register, all
56: .seealso: TSGLAdaptRegisterDestroy()
57: @*/
58: PetscErrorCode TSGLAdaptRegisterAll(const char path[])
59: {
63: TSGLAdaptRegisterDynamic(TSGLADAPT_NONE,path,"TSGLAdaptCreate_None",TSGLAdaptCreate_None);
64: TSGLAdaptRegisterDynamic(TSGLADAPT_SIZE,path,"TSGLAdaptCreate_Size",TSGLAdaptCreate_Size);
65: TSGLAdaptRegisterDynamic(TSGLADAPT_BOTH,path,"TSGLAdaptCreate_Both",TSGLAdaptCreate_Both);
66: return(0);
67: }
71: /*@C
72: TSGLFinalizePackage - This function destroys everything in the TSGL package. It is
73: called from PetscFinalize().
75: Level: developer
77: .keywords: Petsc, destroy, package
78: .seealso: PetscFinalize()
79: @*/
80: PetscErrorCode TSGLAdaptFinalizePackage(void)
81: {
83: TSGLAdaptPackageInitialized = PETSC_FALSE;
84: TSGLAdaptRegisterAllCalled = PETSC_FALSE;
85: TSGLAdaptList = PETSC_NULL;
86: return(0);
87: }
91: /*@C
92: TSGLAdaptInitializePackage - This function initializes everything in the TSGLAdapt package. It is
93: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
94: TSCreate_GL() when using static libraries.
96: Input Parameter:
97: path - The dynamic library path, or PETSC_NULL
99: Level: developer
101: .keywords: TSGLAdapt, initialize, package
102: .seealso: PetscInitialize()
103: @*/
104: PetscErrorCode TSGLAdaptInitializePackage(const char path[])
105: {
109: if (TSGLAdaptPackageInitialized) return(0);
110: TSGLAdaptPackageInitialized = PETSC_TRUE;
111: PetscClassIdRegister("TSGLAdapt",&TSGLADAPT_CLASSID);
112: TSGLAdaptRegisterAll(path);
113: PetscRegisterFinalize(TSGLAdaptFinalizePackage);
114: return(0);
115: }
119: /*@C
120: TSGLAdaptRegisterDestroy - Frees the list of adaptivity schemes that were registered by TSGLAdaptRegister()/TSGLAdaptRegisterDynamic().
122: Not Collective
124: Level: advanced
126: .keywords: TSGLAdapt, register, destroy
127: .seealso: TSGLAdaptRegister(), TSGLAdaptRegisterAll(), TSGLAdaptRegisterDynamic()
128: @*/
129: PetscErrorCode TSGLAdaptRegisterDestroy(void)
130: {
134: PetscFListDestroy(&TSGLAdaptList);
135: TSGLAdaptRegisterAllCalled = PETSC_FALSE;
136: return(0);
137: }
142: PetscErrorCode TSGLAdaptSetType(TSGLAdapt adapt,const TSGLAdaptType type)
143: {
144: PetscErrorCode ierr,(*r)(TSGLAdapt);
147: PetscFListFind(TSGLAdaptList,((PetscObject)adapt)->comm,type,PETSC_TRUE,(void(**)(void))&r);
148: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLAdapt type \"%s\" given",type);
149: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
150: (*r)(adapt);
151: PetscObjectChangeTypeName((PetscObject)adapt,type);
152: return(0);
153: }
157: PetscErrorCode TSGLAdaptSetOptionsPrefix(TSGLAdapt adapt,const char prefix[])
158: {
162: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
163: return(0);
164: }
168: PetscErrorCode TSGLAdaptView(TSGLAdapt adapt,PetscViewer viewer)
169: {
171: PetscBool iascii;
174: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
175: if (iascii) {
176: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer,"TSGLAdapt Object");
177: if (adapt->ops->view) {
178: PetscViewerASCIIPushTab(viewer);
179: (*adapt->ops->view)(adapt,viewer);
180: PetscViewerASCIIPopTab(viewer);
181: }
182: } else {
183: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
184: }
186: return(0);
187: }
191: PetscErrorCode TSGLAdaptDestroy(TSGLAdapt *adapt)
192: {
196: if (!*adapt) return(0);
198: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
199: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
200: PetscHeaderDestroy(adapt);
201: return(0);
202: }
206: PetscErrorCode TSGLAdaptSetFromOptions(TSGLAdapt adapt)
207: {
209: char type[256] = TSGLADAPT_BOTH;
210: PetscBool flg;
213: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGL, but currently this
214: * function can only be called from inside TSSetFromOptions_GL() */
215: PetscOptionsHead("TSGL Adaptivity options");
216: PetscOptionsList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLAdaptSetType",TSGLAdaptList,
217: ((PetscObject)adapt)->type_name?((PetscObject)adapt)->type_name:type,type,sizeof type,&flg);
218: if (flg || !((PetscObject)adapt)->type_name) {
219: TSGLAdaptSetType(adapt,type);
220: }
221: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(adapt);}
222: PetscOptionsTail();
223: return(0);
224: }
228: 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)
229: {
240: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
241: return(0);
242: }
246: PetscErrorCode TSGLAdaptCreate(MPI_Comm comm,TSGLAdapt *inadapt)
247: {
249: TSGLAdapt adapt;
252: *inadapt = 0;
253: PetscHeaderCreate(adapt,_p_TSGLAdapt,struct _TSGLAdaptOps,TSGLADAPT_CLASSID,0,"TSGLAdapt","General Linear adaptivity","TS",comm,TSGLAdaptDestroy,TSGLAdaptView);
254: *inadapt = adapt;
255: return(0);
256: }
259: /*
260: * Implementations
261: */
265: static PetscErrorCode TSGLAdaptDestroy_JustFree(TSGLAdapt adapt)
266: {
270: PetscFree(adapt->data);
271: return(0);
272: }
274: /* -------------------------------- None ----------------------------------- */
275: typedef struct {
276: PetscInt scheme;
277: PetscReal h;
278: } TSGLAdapt_None;
282: 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)
283: {
286: *next_sc = cur;
287: *next_h = h;
288: if (*next_h > tleft) {
289: *finish = PETSC_TRUE;
290: *next_h = tleft;
291: } else {
292: *finish = PETSC_FALSE;
293: }
294: return(0);
295: }
297: EXTERN_C_BEGIN
300: PetscErrorCode TSGLAdaptCreate_None(TSGLAdapt adapt)
301: {
303: TSGLAdapt_None *a;
306: PetscNewLog(adapt,TSGLAdapt_None,&a);
307: adapt->data = (void*)a;
308: adapt->ops->choose = TSGLAdaptChoose_None;
309: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
310: return(0);
311: }
312: EXTERN_C_END
314: /* -------------------------------- Size ----------------------------------- */
315: typedef struct {
316: PetscReal desired_h;
317: } TSGLAdapt_Size;
322: 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)
323: {
324: TSGLAdapt_Size *sz = (TSGLAdapt_Size*)adapt->data;
325: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
328: *next_sc = cur;
329: optimal = pow((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
330: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the the current step size and the
331: * one that would have been taken (without smoothing) on the last step. */
332: last_desired_h = sz->desired_h;
333: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
334: if (last_desired_h > 1e-14) { /* Normally only happens on the first step */
335: *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
336: } else {
337: *next_h = sz->desired_h;
338: }
339: if (*next_h > tleft) {
340: *finish = PETSC_TRUE;
341: *next_h = tleft;
342: } else {
343: *finish = PETSC_FALSE;
344: }
345: return(0);
346: }
348: EXTERN_C_BEGIN
351: PetscErrorCode TSGLAdaptCreate_Size(TSGLAdapt adapt)
352: {
354: TSGLAdapt_Size *a;
357: PetscNewLog(adapt,TSGLAdapt_Size,&a);
358: adapt->data = (void*)a;
359: adapt->ops->choose = TSGLAdaptChoose_Size;
360: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
361: return(0);
362: }
363: EXTERN_C_END
365: /* -------------------------------- Both ----------------------------------- */
366: typedef struct {
367: PetscInt count_at_order;
368: PetscReal desired_h;
369: } TSGLAdapt_Both;
374: 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)
375: {
376: TSGLAdapt_Both *both = (TSGLAdapt_Both*)adapt->data;
378: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
379: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
380: PetscInt i;
383: for (i=0; i<n; i++) {
384: PetscReal optimal;
385: trial.id = i;
386: optimal = pow((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
387: trial.h = h*optimal;
388: trial.eff = trial.h/cost[i];
389: if (trial.eff > best.eff) {PetscMemcpy(&best,&trial,sizeof(trial));}
390: if (i == cur) {PetscMemcpy(¤t,&trial,sizeof(trial));}
391: }
392: /* Only switch orders if the scheme offers significant benefits over the current one.
393: When the scheme is not changing, only change step size if it offers significant benefits. */
394: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
395: PetscReal last_desired_h;
396: *next_sc = current.id;
397: last_desired_h = both->desired_h;
398: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
399: *next_h = (both->count_at_order > 0)
400: ? PetscSqrtReal(last_desired_h * both->desired_h)
401: : both->desired_h;
402: both->count_at_order++;
403: } else {
404: PetscReal rat = cost[best.id]/cost[cur];
405: *next_sc = best.id;
406: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
407: both->count_at_order = 0;
408: both->desired_h = best.h;
409: }
411: if (*next_h > tleft) {
412: *finish = PETSC_TRUE;
413: *next_h = tleft;
414: } else {
415: *finish = PETSC_FALSE;
416: }
417: return(0);
418: }
420: EXTERN_C_BEGIN
423: PetscErrorCode TSGLAdaptCreate_Both(TSGLAdapt adapt)
424: {
426: TSGLAdapt_Both *a;
429: PetscNewLog(adapt,TSGLAdapt_Both,&a);
430: adapt->data = (void*)a;
431: adapt->ops->choose = TSGLAdaptChoose_Both;
432: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
433: return(0);
434: }
435: EXTERN_C_END