Actual source code: tsadapt.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
4: static PetscFList TSAdaptList;
5: static PetscBool TSAdaptPackageInitialized;
6: static PetscBool TSAdaptRegisterAllCalled;
7: static PetscClassId TSADAPT_CLASSID;
9: EXTERN_C_BEGIN
10: PetscErrorCode TSAdaptCreate_Basic(TSAdapt);
11: PetscErrorCode TSAdaptCreate_None(TSAdapt);
12: PetscErrorCode TSAdaptCreate_CFL(TSAdapt);
13: EXTERN_C_END
17: /*@C
18: TSAdaptRegister - see TSAdaptRegisterDynamic()
20: Level: advanced
21: @*/
22: PetscErrorCode TSAdaptRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(TSAdapt))
23: {
25: char fullname[PETSC_MAX_PATH_LEN];
28: PetscFListConcat(path,name,fullname);
29: PetscFListAdd(&TSAdaptList,sname,fullname,(void(*)(void))function);
30: return(0);
31: }
35: /*@C
36: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
38: Not Collective
40: Level: advanced
42: .keywords: TSAdapt, register, all
44: .seealso: TSAdaptRegisterDestroy()
45: @*/
46: PetscErrorCode TSAdaptRegisterAll(const char path[])
47: {
51: TSAdaptRegisterDynamic(TSADAPTBASIC,path,"TSAdaptCreate_Basic",TSAdaptCreate_Basic);
52: TSAdaptRegisterDynamic(TSADAPTNONE, path,"TSAdaptCreate_None", TSAdaptCreate_None);
53: TSAdaptRegisterDynamic(TSADAPTCFL, path,"TSAdaptCreate_CFL", TSAdaptCreate_CFL);
54: return(0);
55: }
59: /*@C
60: TSFinalizePackage - This function destroys everything in the TS package. It is
61: called from PetscFinalize().
63: Level: developer
65: .keywords: Petsc, destroy, package
66: .seealso: PetscFinalize()
67: @*/
68: PetscErrorCode TSAdaptFinalizePackage(void)
69: {
71: TSAdaptPackageInitialized = PETSC_FALSE;
72: TSAdaptRegisterAllCalled = PETSC_FALSE;
73: TSAdaptList = PETSC_NULL;
74: return(0);
75: }
79: /*@C
80: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
81: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
82: TSCreate_GL() when using static libraries.
84: Input Parameter:
85: path - The dynamic library path, or PETSC_NULL
87: Level: developer
89: .keywords: TSAdapt, initialize, package
90: .seealso: PetscInitialize()
91: @*/
92: PetscErrorCode TSAdaptInitializePackage(const char path[])
93: {
97: if (TSAdaptPackageInitialized) return(0);
98: TSAdaptPackageInitialized = PETSC_TRUE;
99: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
100: TSAdaptRegisterAll(path);
101: PetscRegisterFinalize(TSAdaptFinalizePackage);
102: return(0);
103: }
107: /*@C
108: TSAdaptRegisterDestroy - Frees the list of adaptivity schemes that were registered by TSAdaptRegister()/TSAdaptRegisterDynamic().
110: Not Collective
112: Level: advanced
114: .keywords: TSAdapt, register, destroy
115: .seealso: TSAdaptRegister(), TSAdaptRegisterAll(), TSAdaptRegisterDynamic()
116: @*/
117: PetscErrorCode TSAdaptRegisterDestroy(void)
118: {
122: PetscFListDestroy(&TSAdaptList);
123: TSAdaptRegisterAllCalled = PETSC_FALSE;
124: return(0);
125: }
130: PetscErrorCode TSAdaptSetType(TSAdapt adapt,const TSAdaptType type)
131: {
132: PetscErrorCode ierr,(*r)(TSAdapt);
135: PetscFListFind(TSAdaptList,((PetscObject)adapt)->comm,type,PETSC_TRUE,(void(**)(void))&r);
136: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt 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: }
145: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
146: {
150: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
151: return(0);
152: }
156: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
157: {
159: PetscBool iascii;
162: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
163: if (iascii) {
164: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer,"TSAdapt Object");
165: PetscViewerASCIIPrintf(viewer,"number of candidates %D\n",adapt->candidates.n);
166: if (adapt->ops->view) {
167: PetscViewerASCIIPushTab(viewer);
168: (*adapt->ops->view)(adapt,viewer);
169: PetscViewerASCIIPopTab(viewer);
170: }
171: } else {
172: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported",((PetscObject)viewer)->type_name);
173: }
174: return(0);
175: }
179: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
180: {
184: if (!*adapt) return(0);
186: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
187: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
188: PetscViewerDestroy(&(*adapt)->monitor);
189: PetscHeaderDestroy(adapt);
190: return(0);
191: }
195: /*@
196: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
198: Collective on TSAdapt
200: Input Arguments:
201: + adapt - adaptive controller context
202: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
204: Level: intermediate
206: .seealso: TSAdaptChoose()
207: @*/
208: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
209: {
213: if (flg) {
214: if (!adapt->monitor) {PetscViewerASCIIOpen(((PetscObject)adapt)->comm,"stdout",&adapt->monitor);}
215: } else {
216: PetscViewerDestroy(&adapt->monitor);
217: }
218: return(0);
219: }
223: /*@C
224: TSAdaptSetCheckStage - set a callback to check convergence for a stage
226: Logically collective on TSAdapt
228: Input Arguments:
229: + adapt - adaptive controller context
230: - func - stage check function
232: Arguments of func:
233: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
235: + adapt - adaptive controller context
236: . ts - time stepping context
237: - accept - pending choice of whether to accept, can be modified by this routine
239: Level: advanced
241: .seealso: TSAdaptChoose()
242: @*/
243: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscBool*))
244: {
248: adapt->ops->checkstage = func;
249: return(0);
250: }
254: /*@
255: TSAdaptSetStepLimits - Set minimum and maximum step sizes to be considered by the controller
257: Logically Collective
259: Input Arguments:
260: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
261: . hmin - minimum time step
262: - hmax - maximum time step
264: Options Database Keys:
265: + -ts_adapt_dt_min - minimum time step
266: - -ts_adapt_dt_max - maximum time step
268: Level: intermediate
270: .seealso: TSAdapt
271: @*/
272: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
273: {
276: if (hmin != PETSC_DECIDE) adapt->dt_min = hmin;
277: if (hmax != PETSC_DECIDE) adapt->dt_max = hmax;
278: return(0);
279: }
283: /*@
284: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
286: Collective on TSAdapt
288: Input Parameter:
289: . adapt - the TSAdapt context
291: Options Database Keys:
292: . -ts_adapt_type <type> - basic
294: Level: advanced
296: Notes:
297: This function is automatically called by TSSetFromOptions()
299: .keywords: TS, TSGetAdapt(), TSAdaptSetType()
301: .seealso: TSGetType()
302: @*/
303: PetscErrorCode TSAdaptSetFromOptions(TSAdapt adapt)
304: {
306: char type[256] = TSADAPTBASIC;
307: PetscBool set,flg;
310: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
311: * function can only be called from inside TSSetFromOptions_GL() */
312: PetscOptionsHead("TS Adaptivity options");
313: PetscOptionsList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,
314: ((PetscObject)adapt)->type_name?((PetscObject)adapt)->type_name:type,type,sizeof type,&flg);
315: if (flg || !((PetscObject)adapt)->type_name) {
316: TSAdaptSetType(adapt,type);
317: }
318: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(adapt);}
319: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",adapt->dt_min,&adapt->dt_min,PETSC_NULL);
320: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",adapt->dt_max,&adapt->dt_max,PETSC_NULL);
321: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,PETSC_NULL);
322: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
323: if (set) {TSAdaptSetMonitor(adapt,flg);}
324: PetscOptionsTail();
325: return(0);
326: }
330: /*@
331: TSAdaptCandidatesClear - clear any previously set candidate schemes
333: Logically Collective
335: Input Argument:
336: . adapt - adaptive controller
338: Level: developer
340: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
341: @*/
342: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
343: {
347: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
348: return(0);
349: }
353: /*@C
354: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
356: Logically Collective
358: Input Arguments:
359: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
360: . name - name of the candidate scheme to add
361: . order - order of the candidate scheme
362: . stageorder - stage order of the candidate scheme
363: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
364: . cost - relative measure of the amount of work required for the candidate scheme
365: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
367: Note:
368: This routine is not available in Fortran.
370: Level: developer
372: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
373: @*/
374: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
375: {
376: PetscInt c;
380: if (order < 1) SETERRQ1(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
381: if (inuse) {
382: if (adapt->candidates.inuse_set) SETERRQ(((PetscObject)adapt)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
383: adapt->candidates.inuse_set = PETSC_TRUE;
384: }
385: /* first slot if this is the current scheme, otherwise the next available slot */
386: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
387: adapt->candidates.name[c] = name;
388: adapt->candidates.order[c] = order;
389: adapt->candidates.stageorder[c] = stageorder;
390: adapt->candidates.ccfl[c] = ccfl;
391: adapt->candidates.cost[c] = cost;
392: adapt->candidates.n++;
393: return(0);
394: }
398: /*@C
399: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
401: Not Collective
403: Input Arguments:
404: . adapt - time step adaptivity context
406: Output Arguments:
407: + n - number of candidate schemes, always at least 1
408: . order - the order of each candidate scheme
409: . stageorder - the stage order of each candidate scheme
410: . ccfl - the CFL coefficient of each scheme
411: - cost - the relative cost of each scheme
413: Level: developer
415: Note:
416: The current scheme is always returned in the first slot
418: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
419: @*/
420: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
421: {
424: if (n) *n = adapt->candidates.n;
425: if (order) *order = adapt->candidates.order;
426: if (stageorder) *stageorder = adapt->candidates.stageorder;
427: if (ccfl) *ccfl = adapt->candidates.ccfl;
428: if (cost) *cost = adapt->candidates.cost;
429: return(0);
430: }
434: /*@C
435: TSAdaptChoose - choose which method and step size to use for the next step
437: Logically Collective
439: Input Arguments:
440: + adapt - adaptive contoller
441: - h - current step size
443: Output Arguments:
444: + next_sc - scheme to use for the next step
445: . next_h - step size to use for the next step
446: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
448: Note:
449: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
450: being retried after an initial rejection.
452: Level: developer
454: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
455: @*/
456: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
457: {
459: PetscReal wlte = -1.0;
467: if (adapt->candidates.n < 1) SETERRQ1(((PetscObject)adapt)->comm,PETSC_ERR_ARG_WRONGSTATE,"%D candidates have been registered",adapt->candidates.n);
468: if (!adapt->candidates.inuse_set) SETERRQ1(((PetscObject)adapt)->comm,PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n);
469: (*adapt->ops->choose)(adapt,ts,h,next_sc,next_h,accept,&wlte);
470: if (*next_sc < 0 || adapt->candidates.n <= *next_sc) SETERRQ2(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",*next_sc,adapt->candidates.n-1);
471: if (!(*next_h > 0.)) SETERRQ1(((PetscObject)adapt)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %G must be positive",*next_h);
473: if (adapt->monitor) {
474: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
475: if (wlte < 0) {
476: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D %s t=%-11g+%10.3e family='%s' scheme=%D:'%s' dt=%-10g\n",((PetscObject)adapt)->type_name,ts->steps,*accept?"accepted":"rejected",(double)ts->ptime,(double)h,((PetscObject)ts)->type_name,*next_sc,adapt->candidates.name[*next_sc],(double)*next_h);
477: } else {
478: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D %s t=%-11g+%10.3e wlte=%5.3g family='%s' scheme=%D:'%s' dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,*accept?"accepted":"rejected",(double)ts->ptime,(double)h,(double)wlte,((PetscObject)ts)->type_name,*next_sc,adapt->candidates.name[*next_sc],(double)*next_h);
479: }
480: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
481: }
482: return(0);
483: }
487: /*@
488: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)
490: Collective
492: Input Arguments:
493: + adapt - adaptive controller context
494: - ts - time stepper
496: Output Arguments:
497: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
499: Level: developer
501: .seealso:
502: @*/
503: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscBool *accept)
504: {
505: PetscErrorCode ierr;
506: SNES snes;
507: SNESConvergedReason snesreason;
510: *accept = PETSC_TRUE;
511: TSGetSNES(ts,&snes);
512: SNESGetConvergedReason(snes,&snesreason);
513: if (snesreason < 0) {
514: PetscReal dt,new_dt;
515: *accept = PETSC_FALSE;
516: TSGetTimeStep(ts,&dt);
517: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
518: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
519: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
520: if (adapt->monitor) {
521: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
522: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D stage rejected t=%-11g+%10.3e, %D failures exceeds current TS allowed\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,dt,ts->num_snes_failures);
523: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
524: }
525: } else {
526: new_dt = dt*adapt->scale_solve_failed;
527: TSSetTimeStep(ts,new_dt);
528: if (adapt->monitor) {
529: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
530: PetscViewerASCIIPrintf(adapt->monitor," TSAdapt '%s': step %3D stage rejected t=%-11g+%10.3e retrying with dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)dt,(double)new_dt);
531: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
532: }
533: }
534: }
535: if (adapt->ops->checkstage) {(*adapt->ops->checkstage)(adapt,ts,accept);}
536: return(0);
537: }
543: /*@
544: TSAdaptCreate - create an adaptive controller context for time stepping
546: Collective on MPI_Comm
548: Input Parameter:
549: . comm - The communicator
551: Output Parameter:
552: . adapt - new TSAdapt object
554: Level: developer
556: Notes:
557: TSAdapt creation is handled by TS, so users should not need to call this function.
559: .keywords: TSAdapt, create
560: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
561: @*/
562: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
563: {
565: TSAdapt adapt;
568: *inadapt = 0;
569: PetscHeaderCreate(adapt,_p_TSAdapt,struct _TSAdaptOps,TSADAPT_CLASSID,0,"TSAdapt","General Linear adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
571: adapt->dt_min = 1e-20;
572: adapt->dt_max = 1e50;
573: adapt->scale_solve_failed = 0.25;
575: *inadapt = adapt;
576: return(0);
577: }