Actual source code: gladapt.c
petsc-3.6.1 2015-08-06
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)(PetscOptions*,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: if (TSGLAdaptRegisterAllCalled) return(0);
83: TSGLAdaptRegisterAllCalled = PETSC_TRUE;
84: TSGLAdaptRegister(TSGLADAPT_NONE,TSGLAdaptCreate_None);
85: TSGLAdaptRegister(TSGLADAPT_SIZE,TSGLAdaptCreate_Size);
86: TSGLAdaptRegister(TSGLADAPT_BOTH,TSGLAdaptCreate_Both);
87: return(0);
88: }
92: /*@C
93: TSGLFinalizePackage - This function destroys everything in the TSGL package. It is
94: called from PetscFinalize().
96: Level: developer
98: .keywords: Petsc, destroy, package
99: .seealso: PetscFinalize()
100: @*/
101: PetscErrorCode TSGLAdaptFinalizePackage(void)
102: {
106: PetscFunctionListDestroy(&TSGLAdaptList);
107: TSGLAdaptPackageInitialized = PETSC_FALSE;
108: TSGLAdaptRegisterAllCalled = PETSC_FALSE;
109: return(0);
110: }
114: /*@C
115: TSGLAdaptInitializePackage - This function initializes everything in the TSGLAdapt package. It is
116: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
117: TSCreate_GL() when using static libraries.
119: Level: developer
121: .keywords: TSGLAdapt, initialize, package
122: .seealso: PetscInitialize()
123: @*/
124: PetscErrorCode TSGLAdaptInitializePackage(void)
125: {
129: if (TSGLAdaptPackageInitialized) return(0);
130: TSGLAdaptPackageInitialized = PETSC_TRUE;
131: PetscClassIdRegister("TSGLAdapt",&TSGLADAPT_CLASSID);
132: TSGLAdaptRegisterAll();
133: PetscRegisterFinalize(TSGLAdaptFinalizePackage);
134: return(0);
135: }
139: PetscErrorCode TSGLAdaptSetType(TSGLAdapt adapt,TSGLAdaptType type)
140: {
141: PetscErrorCode ierr,(*r)(TSGLAdapt);
144: PetscFunctionListFind(TSGLAdaptList,type,&r);
145: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLAdapt type \"%s\" given",type);
146: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
147: (*r)(adapt);
148: PetscObjectChangeTypeName((PetscObject)adapt,type);
149: return(0);
150: }
154: PetscErrorCode TSGLAdaptSetOptionsPrefix(TSGLAdapt adapt,const char prefix[])
155: {
159: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
160: return(0);
161: }
165: PetscErrorCode TSGLAdaptView(TSGLAdapt adapt,PetscViewer viewer)
166: {
168: PetscBool iascii;
171: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
172: if (iascii) {
173: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
174: if (adapt->ops->view) {
175: PetscViewerASCIIPushTab(viewer);
176: (*adapt->ops->view)(adapt,viewer);
177: PetscViewerASCIIPopTab(viewer);
178: }
179: }
180: return(0);
181: }
185: PetscErrorCode TSGLAdaptDestroy(TSGLAdapt *adapt)
186: {
190: if (!*adapt) return(0);
192: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
193: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
194: PetscHeaderDestroy(adapt);
195: return(0);
196: }
200: PetscErrorCode TSGLAdaptSetFromOptions(PetscOptions *PetscOptionsObject,TSGLAdapt adapt)
201: {
203: char type[256] = TSGLADAPT_BOTH;
204: PetscBool flg;
207: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGL, but currently this
208: * function can only be called from inside TSSetFromOptions_GL() */
209: PetscOptionsHead(PetscOptionsObject,"TSGL Adaptivity options");
210: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLAdaptSetType",TSGLAdaptList,
211: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
212: if (flg || !((PetscObject)adapt)->type_name) {
213: TSGLAdaptSetType(adapt,type);
214: }
215: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
216: PetscOptionsTail();
217: return(0);
218: }
222: 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)
223: {
234: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
235: return(0);
236: }
240: PetscErrorCode TSGLAdaptCreate(MPI_Comm comm,TSGLAdapt *inadapt)
241: {
243: TSGLAdapt adapt;
246: *inadapt = NULL;
247: PetscHeaderCreate(adapt,TSGLADAPT_CLASSID,"TSGLAdapt","General Linear adaptivity","TS",comm,TSGLAdaptDestroy,TSGLAdaptView);
248: *inadapt = adapt;
249: return(0);
250: }
253: /*
254: * Implementations
255: */
259: static PetscErrorCode TSGLAdaptDestroy_JustFree(TSGLAdapt adapt)
260: {
264: PetscFree(adapt->data);
265: return(0);
266: }
268: /* -------------------------------- None ----------------------------------- */
269: typedef struct {
270: PetscInt scheme;
271: PetscReal h;
272: } TSGLAdapt_None;
276: 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)
277: {
280: *next_sc = cur;
281: *next_h = h;
282: if (*next_h > tleft) {
283: *finish = PETSC_TRUE;
284: *next_h = tleft;
285: } else *finish = PETSC_FALSE;
286: return(0);
287: }
291: PetscErrorCode TSGLAdaptCreate_None(TSGLAdapt adapt)
292: {
294: TSGLAdapt_None *a;
297: PetscNewLog(adapt,&a);
298: adapt->data = (void*)a;
299: adapt->ops->choose = TSGLAdaptChoose_None;
300: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
301: return(0);
302: }
304: /* -------------------------------- Size ----------------------------------- */
305: typedef struct {
306: PetscReal desired_h;
307: } TSGLAdapt_Size;
312: 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)
313: {
314: TSGLAdapt_Size *sz = (TSGLAdapt_Size*)adapt->data;
315: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
318: *next_sc = cur;
319: optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
320: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
321: * one that would have been taken (without smoothing) on the last step. */
322: last_desired_h = sz->desired_h;
323: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
325: /* Normally only happens on the first step */
326: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
327: else *next_h = sz->desired_h;
329: if (*next_h > tleft) {
330: *finish = PETSC_TRUE;
331: *next_h = tleft;
332: } else *finish = PETSC_FALSE;
333: return(0);
334: }
338: PetscErrorCode TSGLAdaptCreate_Size(TSGLAdapt adapt)
339: {
341: TSGLAdapt_Size *a;
344: PetscNewLog(adapt,&a);
345: adapt->data = (void*)a;
346: adapt->ops->choose = TSGLAdaptChoose_Size;
347: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
348: return(0);
349: }
351: /* -------------------------------- Both ----------------------------------- */
352: typedef struct {
353: PetscInt count_at_order;
354: PetscReal desired_h;
355: } TSGLAdapt_Both;
360: 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)
361: {
362: TSGLAdapt_Both *both = (TSGLAdapt_Both*)adapt->data;
364: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
365: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
366: PetscInt i;
369: for (i=0; i<n; i++) {
370: PetscReal optimal;
371: trial.id = i;
372: optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
373: trial.h = h*optimal;
374: trial.eff = trial.h/cost[i];
375: if (trial.eff > best.eff) {PetscMemcpy(&best,&trial,sizeof(trial));}
376: if (i == cur) {PetscMemcpy(¤t,&trial,sizeof(trial));}
377: }
378: /* Only switch orders if the scheme offers significant benefits over the current one.
379: When the scheme is not changing, only change step size if it offers significant benefits. */
380: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
381: PetscReal last_desired_h;
382: *next_sc = current.id;
383: last_desired_h = both->desired_h;
384: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
385: *next_h = (both->count_at_order > 0)
386: ? PetscSqrtReal(last_desired_h * both->desired_h)
387: : both->desired_h;
388: both->count_at_order++;
389: } else {
390: PetscReal rat = cost[best.id]/cost[cur];
391: *next_sc = best.id;
392: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
393: both->count_at_order = 0;
394: both->desired_h = best.h;
395: }
397: if (*next_h > tleft) {
398: *finish = PETSC_TRUE;
399: *next_h = tleft;
400: } else *finish = PETSC_FALSE;
401: return(0);
402: }
406: PetscErrorCode TSGLAdaptCreate_Both(TSGLAdapt adapt)
407: {
409: TSGLAdapt_Both *a;
412: PetscNewLog(adapt,&a);
413: adapt->data = (void*)a;
414: adapt->ops->choose = TSGLAdaptChoose_Both;
415: adapt->ops->destroy = TSGLAdaptDestroy_JustFree;
416: return(0);
417: }