Actual source code: glleadapt.c
petsc-3.10.5 2019-03-28
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: .keywords: TSGLLEAdapt, register
51: .seealso: TSGLLEAdaptRegisterAll()
52: @*/
53: PetscErrorCode TSGLLEAdaptRegister(const char sname[],PetscErrorCode (*function)(TSGLLEAdapt))
54: {
58: TSGLLEAdaptInitializePackage();
59: PetscFunctionListAdd(&TSGLLEAdaptList,sname,function);
60: return(0);
61: }
63: /*@C
64: TSGLLEAdaptRegisterAll - Registers all of the adaptivity schemes in TSGLLEAdapt
66: Not Collective
68: Level: advanced
70: .keywords: TSGLLEAdapt, register, all
72: .seealso: TSGLLEAdaptRegisterDestroy()
73: @*/
74: PetscErrorCode TSGLLEAdaptRegisterAll(void)
75: {
79: if (TSGLLEAdaptRegisterAllCalled) return(0);
80: TSGLLEAdaptRegisterAllCalled = PETSC_TRUE;
81: TSGLLEAdaptRegister(TSGLLEADAPT_NONE,TSGLLEAdaptCreate_None);
82: TSGLLEAdaptRegister(TSGLLEADAPT_SIZE,TSGLLEAdaptCreate_Size);
83: TSGLLEAdaptRegister(TSGLLEADAPT_BOTH,TSGLLEAdaptCreate_Both);
84: return(0);
85: }
87: /*@C
88: TSGLLEFinalizePackage - This function destroys everything in the TSGLLE package. It is
89: called from PetscFinalize().
91: Level: developer
93: .keywords: Petsc, destroy, package
94: .seealso: PetscFinalize()
95: @*/
96: PetscErrorCode TSGLLEAdaptFinalizePackage(void)
97: {
101: PetscFunctionListDestroy(&TSGLLEAdaptList);
102: TSGLLEAdaptPackageInitialized = PETSC_FALSE;
103: TSGLLEAdaptRegisterAllCalled = PETSC_FALSE;
104: return(0);
105: }
107: /*@C
108: TSGLLEAdaptInitializePackage - This function initializes everything in the TSGLLEAdapt package. It is
109: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
110: TSCreate_GLLE() when using static libraries.
112: Level: developer
114: .keywords: TSGLLEAdapt, initialize, package
115: .seealso: PetscInitialize()
116: @*/
117: PetscErrorCode TSGLLEAdaptInitializePackage(void)
118: {
122: if (TSGLLEAdaptPackageInitialized) return(0);
123: TSGLLEAdaptPackageInitialized = PETSC_TRUE;
124: PetscClassIdRegister("TSGLLEAdapt",&TSGLLEADAPT_CLASSID);
125: TSGLLEAdaptRegisterAll();
126: PetscRegisterFinalize(TSGLLEAdaptFinalizePackage);
127: return(0);
128: }
130: PetscErrorCode TSGLLEAdaptSetType(TSGLLEAdapt adapt,TSGLLEAdaptType type)
131: {
132: PetscErrorCode ierr,(*r)(TSGLLEAdapt);
135: PetscFunctionListFind(TSGLLEAdaptList,type,&r);
136: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSGLLEAdapt type \"%s\" given",type);
137: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
138: (*r)(adapt);
139: PetscObjectChangeTypeName((PetscObject)adapt,type);
140: return(0);
141: }
143: PetscErrorCode TSGLLEAdaptSetOptionsPrefix(TSGLLEAdapt adapt,const char prefix[])
144: {
148: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
149: return(0);
150: }
152: PetscErrorCode TSGLLEAdaptView(TSGLLEAdapt adapt,PetscViewer viewer)
153: {
155: PetscBool iascii;
158: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
159: if (iascii) {
160: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);
161: if (adapt->ops->view) {
162: PetscViewerASCIIPushTab(viewer);
163: (*adapt->ops->view)(adapt,viewer);
164: PetscViewerASCIIPopTab(viewer);
165: }
166: }
167: return(0);
168: }
170: PetscErrorCode TSGLLEAdaptDestroy(TSGLLEAdapt *adapt)
171: {
175: if (!*adapt) return(0);
177: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
178: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
179: PetscHeaderDestroy(adapt);
180: return(0);
181: }
183: PetscErrorCode TSGLLEAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSGLLEAdapt adapt)
184: {
186: char type[256] = TSGLLEADAPT_BOTH;
187: PetscBool flg;
190: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TSGLLE, but currently this
191: * function can only be called from inside TSSetFromOptions_GLLE() */
192: PetscOptionsHead(PetscOptionsObject,"TSGLLE Adaptivity options");
193: PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSGLLEAdaptSetType",TSGLLEAdaptList,
194: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
195: if (flg || !((PetscObject)adapt)->type_name) {
196: TSGLLEAdaptSetType(adapt,type);
197: }
198: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);}
199: PetscOptionsTail();
200: return(0);
201: }
203: 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)
204: {
215: (*adapt->ops->choose)(adapt,n,orders,errors,cost,cur,h,tleft,next_sc,next_h,finish);
216: return(0);
217: }
219: PetscErrorCode TSGLLEAdaptCreate(MPI_Comm comm,TSGLLEAdapt *inadapt)
220: {
222: TSGLLEAdapt adapt;
225: *inadapt = NULL;
226: PetscHeaderCreate(adapt,TSGLLEADAPT_CLASSID,"TSGLLEAdapt","General Linear adaptivity","TS",comm,TSGLLEAdaptDestroy,TSGLLEAdaptView);
227: *inadapt = adapt;
228: return(0);
229: }
232: /*
233: * Implementations
234: */
236: static PetscErrorCode TSGLLEAdaptDestroy_JustFree(TSGLLEAdapt adapt)
237: {
241: PetscFree(adapt->data);
242: return(0);
243: }
245: /* -------------------------------- None ----------------------------------- */
246: typedef struct {
247: PetscInt scheme;
248: PetscReal h;
249: } TSGLLEAdapt_None;
251: 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)
252: {
255: *next_sc = cur;
256: *next_h = h;
257: if (*next_h > tleft) {
258: *finish = PETSC_TRUE;
259: *next_h = tleft;
260: } else *finish = PETSC_FALSE;
261: return(0);
262: }
264: PetscErrorCode TSGLLEAdaptCreate_None(TSGLLEAdapt adapt)
265: {
267: TSGLLEAdapt_None *a;
270: PetscNewLog(adapt,&a);
271: adapt->data = (void*)a;
272: adapt->ops->choose = TSGLLEAdaptChoose_None;
273: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
274: return(0);
275: }
277: /* -------------------------------- Size ----------------------------------- */
278: typedef struct {
279: PetscReal desired_h;
280: } TSGLLEAdapt_Size;
283: 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)
284: {
285: TSGLLEAdapt_Size *sz = (TSGLLEAdapt_Size*)adapt->data;
286: PetscReal dec = 0.2,inc = 5.0,safe = 0.9,optimal,last_desired_h;
289: *next_sc = cur;
290: optimal = PetscPowReal((PetscReal)errors[cur],(PetscReal)-1./(safe*orders[cur]));
291: /* Step sizes oscillate when there is no smoothing. Here we use a geometric mean of the current step size and the
292: * one that would have been taken (without smoothing) on the last step. */
293: last_desired_h = sz->desired_h;
294: sz->desired_h = h*PetscMax(dec,PetscMin(inc,optimal)); /* Trim to [dec,inc] */
296: /* Normally only happens on the first step */
297: if (last_desired_h > 1e-14) *next_h = PetscSqrtReal(last_desired_h * sz->desired_h);
298: else *next_h = sz->desired_h;
300: if (*next_h > tleft) {
301: *finish = PETSC_TRUE;
302: *next_h = tleft;
303: } else *finish = PETSC_FALSE;
304: return(0);
305: }
307: PetscErrorCode TSGLLEAdaptCreate_Size(TSGLLEAdapt adapt)
308: {
310: TSGLLEAdapt_Size *a;
313: PetscNewLog(adapt,&a);
314: adapt->data = (void*)a;
315: adapt->ops->choose = TSGLLEAdaptChoose_Size;
316: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
317: return(0);
318: }
320: /* -------------------------------- Both ----------------------------------- */
321: typedef struct {
322: PetscInt count_at_order;
323: PetscReal desired_h;
324: } TSGLLEAdapt_Both;
327: 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)
328: {
329: TSGLLEAdapt_Both *both = (TSGLLEAdapt_Both*)adapt->data;
331: PetscReal dec = 0.2,inc = 5.0,safe = 0.9;
332: struct {PetscInt id; PetscReal h,eff;} best={-1,0,0},trial={-1,0,0},current={-1,0,0};
333: PetscInt i;
336: for (i=0; i<n; i++) {
337: PetscReal optimal;
338: trial.id = i;
339: optimal = PetscPowReal((PetscReal)errors[i],(PetscReal)-1./(safe*orders[i]));
340: trial.h = h*optimal;
341: trial.eff = trial.h/cost[i];
342: if (trial.eff > best.eff) {PetscMemcpy(&best,&trial,sizeof(trial));}
343: if (i == cur) {PetscMemcpy(¤t,&trial,sizeof(trial));}
344: }
345: /* Only switch orders if the scheme offers significant benefits over the current one.
346: When the scheme is not changing, only change step size if it offers significant benefits. */
347: if (best.eff < 1.2*current.eff || both->count_at_order < orders[cur]+2) {
348: PetscReal last_desired_h;
349: *next_sc = current.id;
350: last_desired_h = both->desired_h;
351: both->desired_h = PetscMax(h*dec,PetscMin(h*inc,current.h));
352: *next_h = (both->count_at_order > 0)
353: ? PetscSqrtReal(last_desired_h * both->desired_h)
354: : both->desired_h;
355: both->count_at_order++;
356: } else {
357: PetscReal rat = cost[best.id]/cost[cur];
358: *next_sc = best.id;
359: *next_h = PetscMax(h*rat*dec,PetscMin(h*rat*inc,best.h));
360: both->count_at_order = 0;
361: both->desired_h = best.h;
362: }
364: if (*next_h > tleft) {
365: *finish = PETSC_TRUE;
366: *next_h = tleft;
367: } else *finish = PETSC_FALSE;
368: return(0);
369: }
371: PetscErrorCode TSGLLEAdaptCreate_Both(TSGLLEAdapt adapt)
372: {
374: TSGLLEAdapt_Both *a;
377: PetscNewLog(adapt,&a);
378: adapt->data = (void*)a;
379: adapt->ops->choose = TSGLLEAdaptChoose_Both;
380: adapt->ops->destroy = TSGLLEAdaptDestroy_JustFree;
381: return(0);
382: }