Actual source code: tsadapt.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/tsimpl.h> /*I "petscts.h" I*/
4: static PetscFunctionList TSAdaptList;
5: static PetscBool TSAdaptPackageInitialized;
6: static PetscBool TSAdaptRegisterAllCalled;
7: static PetscClassId TSADAPT_CLASSID;
9: PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt);
10: PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt);
11: PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt);
15: /*@C
16: TSAdaptRegister - adds a TSAdapt implementation
18: Not Collective
20: Input Parameters:
21: + name_scheme - name of user-defined adaptivity scheme
22: - routine_create - routine to create method context
24: Notes:
25: TSAdaptRegister() may be called multiple times to add several user-defined families.
27: Sample usage:
28: .vb
29: TSAdaptRegister("my_scheme",MySchemeCreate);
30: .ve
32: Then, your scheme can be chosen with the procedural interface via
33: $ TSAdaptSetType(ts,"my_scheme")
34: or at runtime via the option
35: $ -ts_adapt_type my_scheme
37: Level: advanced
39: .keywords: TSAdapt, register
41: .seealso: TSAdaptRegisterAll()
42: @*/
43: PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt))
44: {
48: PetscFunctionListAdd(&TSAdaptList,sname,function);
49: return(0);
50: }
54: /*@C
55: TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt
57: Not Collective
59: Level: advanced
61: .keywords: TSAdapt, register, all
63: .seealso: TSAdaptRegisterDestroy()
64: @*/
65: PetscErrorCode TSAdaptRegisterAll(void)
66: {
70: TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);
71: TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);
72: TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);
73: return(0);
74: }
78: /*@C
79: TSFinalizePackage - This function destroys everything in the TS package. It is
80: called from PetscFinalize().
82: Level: developer
84: .keywords: Petsc, destroy, package
85: .seealso: PetscFinalize()
86: @*/
87: PetscErrorCode TSAdaptFinalizePackage(void)
88: {
92: PetscFunctionListDestroy(&TSAdaptList);
93: TSAdaptPackageInitialized = PETSC_FALSE;
94: TSAdaptRegisterAllCalled = PETSC_FALSE;
95: return(0);
96: }
100: /*@C
101: TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is
102: called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to
103: TSCreate_GL() when using static libraries.
105: Level: developer
107: .keywords: TSAdapt, initialize, package
108: .seealso: PetscInitialize()
109: @*/
110: PetscErrorCode TSAdaptInitializePackage(void)
111: {
115: if (TSAdaptPackageInitialized) return(0);
116: TSAdaptPackageInitialized = PETSC_TRUE;
117: PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);
118: TSAdaptRegisterAll();
119: PetscRegisterFinalize(TSAdaptFinalizePackage);
120: return(0);
121: }
125: PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type)
126: {
127: PetscErrorCode ierr,(*r)(TSAdapt);
130: PetscFunctionListFind(TSAdaptList,type,&r);
131: if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type);
132: if (((PetscObject)adapt)->type_name) {(*adapt->ops->destroy)(adapt);}
133: (*r)(adapt);
134: PetscObjectChangeTypeName((PetscObject)adapt,type);
135: return(0);
136: }
140: PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[])
141: {
145: PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);
146: return(0);
147: }
151: /*@C
152: TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView().
154: Collective on PetscViewer
156: Input Parameters:
157: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or
158: some related function before a call to TSAdaptLoad().
159: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
160: HDF5 file viewer, obtained from PetscViewerHDF5Open()
162: Level: intermediate
164: Notes:
165: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored.
167: Notes for advanced users:
168: Most users should not need to know the details of the binary storage
169: format, since TSAdaptLoad() and TSAdaptView() completely hide these details.
170: But for anyone who's interested, the standard binary matrix storage
171: format is
172: .vb
173: has not yet been determined
174: .ve
176: .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad()
177: @*/
178: PetscErrorCode TSAdaptLoad(TSAdapt tsadapt, PetscViewer viewer)
179: {
181: PetscBool isbinary;
182: char type[256];
187: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
188: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
190: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
191: TSAdaptSetType(tsadapt, type);
192: if (tsadapt->ops->load) {
193: (*tsadapt->ops->load)(tsadapt,viewer);
194: }
195: return(0);
196: }
200: PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer)
201: {
203: PetscBool iascii,isbinary;
206: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
207: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
208: if (iascii) {
209: PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer,"TSAdapt Object");
210: PetscViewerASCIIPrintf(viewer," number of candidates %D\n",adapt->candidates.n);
211: if (adapt->ops->view) {
212: PetscViewerASCIIPushTab(viewer);
213: (*adapt->ops->view)(adapt,viewer);
214: PetscViewerASCIIPopTab(viewer);
215: }
216: } else if (isbinary) {
217: char type[256];
219: /* need to save FILE_CLASS_ID for adapt class */
220: PetscStrncpy(type,((PetscObject)adapt)->type_name,256);
221: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
222: } else if (adapt->ops->view) {
223: (*adapt->ops->view)(adapt,viewer);
224: }
225: return(0);
226: }
230: PetscErrorCode TSAdaptDestroy(TSAdapt *adapt)
231: {
235: if (!*adapt) return(0);
237: if (--((PetscObject)(*adapt))->refct > 0) {*adapt = 0; return(0);}
238: if ((*adapt)->ops->destroy) {(*(*adapt)->ops->destroy)(*adapt);}
239: PetscViewerDestroy(&(*adapt)->monitor);
240: PetscHeaderDestroy(adapt);
241: return(0);
242: }
246: /*@
247: TSAdaptSetMonitor - Monitor the choices made by the adaptive controller
249: Collective on TSAdapt
251: Input Arguments:
252: + adapt - adaptive controller context
253: - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable
255: Level: intermediate
257: .seealso: TSAdaptChoose()
258: @*/
259: PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg)
260: {
264: if (flg) {
265: if (!adapt->monitor) {PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);}
266: } else {
267: PetscViewerDestroy(&adapt->monitor);
268: }
269: return(0);
270: }
274: /*@C
275: TSAdaptSetCheckStage - set a callback to check convergence for a stage
277: Logically collective on TSAdapt
279: Input Arguments:
280: + adapt - adaptive controller context
281: - func - stage check function
283: Arguments of func:
284: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept)
286: + adapt - adaptive controller context
287: . ts - time stepping context
288: - accept - pending choice of whether to accept, can be modified by this routine
290: Level: advanced
292: .seealso: TSAdaptChoose()
293: @*/
294: PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscBool*))
295: {
299: adapt->ops->checkstage = func;
300: return(0);
301: }
305: /*@
306: TSAdaptSetStepLimits - Set minimum and maximum step sizes to be considered by the controller
308: Logically Collective
310: Input Arguments:
311: + adapt - time step adaptivity context, usually gotten with TSGetAdapt()
312: . hmin - minimum time step
313: - hmax - maximum time step
315: Options Database Keys:
316: + -ts_adapt_dt_min - minimum time step
317: - -ts_adapt_dt_max - maximum time step
319: Level: intermediate
321: .seealso: TSAdapt
322: @*/
323: PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax)
324: {
327: if (hmin != PETSC_DECIDE) adapt->dt_min = hmin;
328: if (hmax != PETSC_DECIDE) adapt->dt_max = hmax;
329: return(0);
330: }
334: /*@
335: TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options.
337: Collective on TSAdapt
339: Input Parameter:
340: . adapt - the TSAdapt context
342: Options Database Keys:
343: . -ts_adapt_type <type> - basic
345: Level: advanced
347: Notes:
348: This function is automatically called by TSSetFromOptions()
350: .keywords: TS, TSGetAdapt(), TSAdaptSetType()
352: .seealso: TSGetType()
353: @*/
354: PetscErrorCode TSAdaptSetFromOptions(TSAdapt adapt)
355: {
357: char type[256] = TSADAPTBASIC;
358: PetscBool set,flg;
361: /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this
362: * function can only be called from inside TSSetFromOptions_GL() */
363: PetscOptionsHead("TS Adaptivity options");
364: PetscOptionsList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,
365: ((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);
366: if (flg || !((PetscObject)adapt)->type_name) {
367: TSAdaptSetType(adapt,type);
368: }
369: if (adapt->ops->setfromoptions) {(*adapt->ops->setfromoptions)(adapt);}
370: PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",adapt->dt_min,&adapt->dt_min,NULL);
371: PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",adapt->dt_max,&adapt->dt_max,NULL);
372: PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);
373: PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);
374: if (set) {TSAdaptSetMonitor(adapt,flg);}
375: PetscOptionsTail();
376: return(0);
377: }
381: /*@
382: TSAdaptCandidatesClear - clear any previously set candidate schemes
384: Logically Collective
386: Input Argument:
387: . adapt - adaptive controller
389: Level: developer
391: .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose()
392: @*/
393: PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt)
394: {
398: PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));
399: return(0);
400: }
404: /*@C
405: TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from
407: Logically Collective
409: Input Arguments:
410: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate()
411: . name - name of the candidate scheme to add
412: . order - order of the candidate scheme
413: . stageorder - stage order of the candidate scheme
414: . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints
415: . cost - relative measure of the amount of work required for the candidate scheme
416: - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme
418: Note:
419: This routine is not available in Fortran.
421: Level: developer
423: .seealso: TSAdaptCandidatesClear(), TSAdaptChoose()
424: @*/
425: PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse)
426: {
427: PetscInt c;
431: if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order);
432: if (inuse) {
433: if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()");
434: adapt->candidates.inuse_set = PETSC_TRUE;
435: }
436: /* first slot if this is the current scheme, otherwise the next available slot */
437: c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n;
439: adapt->candidates.name[c] = name;
440: adapt->candidates.order[c] = order;
441: adapt->candidates.stageorder[c] = stageorder;
442: adapt->candidates.ccfl[c] = ccfl;
443: adapt->candidates.cost[c] = cost;
444: adapt->candidates.n++;
445: return(0);
446: }
450: /*@C
451: TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost
453: Not Collective
455: Input Arguments:
456: . adapt - time step adaptivity context
458: Output Arguments:
459: + n - number of candidate schemes, always at least 1
460: . order - the order of each candidate scheme
461: . stageorder - the stage order of each candidate scheme
462: . ccfl - the CFL coefficient of each scheme
463: - cost - the relative cost of each scheme
465: Level: developer
467: Note:
468: The current scheme is always returned in the first slot
470: .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose()
471: @*/
472: PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost)
473: {
476: if (n) *n = adapt->candidates.n;
477: if (order) *order = adapt->candidates.order;
478: if (stageorder) *stageorder = adapt->candidates.stageorder;
479: if (ccfl) *ccfl = adapt->candidates.ccfl;
480: if (cost) *cost = adapt->candidates.cost;
481: return(0);
482: }
486: /*@C
487: TSAdaptChoose - choose which method and step size to use for the next step
489: Logically Collective
491: Input Arguments:
492: + adapt - adaptive contoller
493: - h - current step size
495: Output Arguments:
496: + next_sc - scheme to use for the next step
497: . next_h - step size to use for the next step
498: - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size
500: Note:
501: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is
502: being retried after an initial rejection.
504: Level: developer
506: .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd()
507: @*/
508: PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept)
509: {
511: PetscReal wlte = -1.0;
519: if (adapt->candidates.n < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"%D candidates have been registered",adapt->candidates.n);
520: if (!adapt->candidates.inuse_set) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"The current in-use scheme is not among the %D candidates",adapt->candidates.n);
521: (*adapt->ops->choose)(adapt,ts,h,next_sc,next_h,accept,&wlte);
522: if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) {
523: /* Reduce time step if it overshoots max time */
524: PetscReal max_time = ts->max_time;
525: PetscReal next_dt = 0.0;
526: if (ts->ptime + ts->time_step + *next_h >= max_time) {
527: next_dt = max_time - (ts->ptime + ts->time_step);
528: if (next_dt > PETSC_SMALL) *next_h = next_dt;
529: else ts->reason = TS_CONVERGED_TIME;
530: }
531: }
532: if (*next_sc < 0 || adapt->candidates.n <= *next_sc) SETERRQ2(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",*next_sc,adapt->candidates.n-1);
533: if (!(*next_h > 0.)) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %G must be positive",*next_h);
535: if (adapt->monitor) {
536: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
537: if (wlte < 0) {
538: 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);
539: } else {
540: 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);
541: }
542: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
543: }
544: return(0);
545: }
549: /*@
550: TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails)
552: Collective
554: Input Arguments:
555: + adapt - adaptive controller context
556: - ts - time stepper
558: Output Arguments:
559: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject
561: Level: developer
563: .seealso:
564: @*/
565: PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscBool *accept)
566: {
567: PetscErrorCode ierr;
568: SNES snes;
569: SNESConvergedReason snesreason;
572: *accept = PETSC_TRUE;
573: TSGetSNES(ts,&snes);
574: SNESGetConvergedReason(snes,&snesreason);
575: if (snesreason < 0) {
576: PetscReal dt,new_dt;
577: *accept = PETSC_FALSE;
578: TSGetTimeStep(ts,&dt);
579: if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) {
580: ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
581: PetscInfo2(ts,"Step=%D, nonlinear solve solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);
582: if (adapt->monitor) {
583: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
584: 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);
585: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
586: }
587: } else {
588: new_dt = dt*adapt->scale_solve_failed;
589: TSSetTimeStep(ts,new_dt);
590: if (adapt->monitor) {
591: PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);
592: 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);
593: PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);
594: }
595: }
596: }
597: if (adapt->ops->checkstage) {(*adapt->ops->checkstage)(adapt,ts,accept);}
598: return(0);
599: }
605: /*@
606: TSAdaptCreate - create an adaptive controller context for time stepping
608: Collective on MPI_Comm
610: Input Parameter:
611: . comm - The communicator
613: Output Parameter:
614: . adapt - new TSAdapt object
616: Level: developer
618: Notes:
619: TSAdapt creation is handled by TS, so users should not need to call this function.
621: .keywords: TSAdapt, create
622: .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy()
623: @*/
624: PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt)
625: {
627: TSAdapt adapt;
630: *inadapt = 0;
631: PetscHeaderCreate(adapt,_p_TSAdapt,struct _TSAdaptOps,TSADAPT_CLASSID,"TSAdapt","General Linear adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);
633: adapt->dt_min = 1e-20;
634: adapt->dt_max = 1e50;
635: adapt->scale_solve_failed = 0.25;
637: *inadapt = adapt;
638: return(0);
639: }