#include /*I "petscts.h" I*/ PetscClassId TSADAPT_CLASSID; static PetscFunctionList TSAdaptList; static PetscBool TSAdaptPackageInitialized; static PetscBool TSAdaptRegisterAllCalled; PETSC_EXTERN PetscErrorCode TSAdaptCreate_None(TSAdapt); PETSC_EXTERN PetscErrorCode TSAdaptCreate_Basic(TSAdapt); PETSC_EXTERN PetscErrorCode TSAdaptCreate_DSP(TSAdapt); PETSC_EXTERN PetscErrorCode TSAdaptCreate_CFL(TSAdapt); PETSC_EXTERN PetscErrorCode TSAdaptCreate_GLEE(TSAdapt); /*@C TSAdaptRegister - adds a TSAdapt implementation Not Collective Input Parameters: + name_scheme - name of user-defined adaptivity scheme - routine_create - routine to create method context Notes: TSAdaptRegister() may be called multiple times to add several user-defined families. Sample usage: .vb TSAdaptRegister("my_scheme",MySchemeCreate); .ve Then, your scheme can be chosen with the procedural interface via $ TSAdaptSetType(ts,"my_scheme") or at runtime via the option $ -ts_adapt_type my_scheme Level: advanced .keywords: TSAdapt, register .seealso: TSAdaptRegisterAll() @*/ PetscErrorCode TSAdaptRegister(const char sname[],PetscErrorCode (*function)(TSAdapt)) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFunctionListAdd(&TSAdaptList,sname,function);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C TSAdaptRegisterAll - Registers all of the adaptivity schemes in TSAdapt Not Collective Level: advanced .keywords: TSAdapt, register, all .seealso: TSAdaptRegisterDestroy() @*/ PetscErrorCode TSAdaptRegisterAll(void) { PetscErrorCode ierr; PetscFunctionBegin; if (TSAdaptRegisterAllCalled) PetscFunctionReturn(0); TSAdaptRegisterAllCalled = PETSC_TRUE; ierr = TSAdaptRegister(TSADAPTNONE, TSAdaptCreate_None);CHKERRQ(ierr); ierr = TSAdaptRegister(TSADAPTBASIC,TSAdaptCreate_Basic);CHKERRQ(ierr); ierr = TSAdaptRegister(TSADAPTDSP, TSAdaptCreate_DSP);CHKERRQ(ierr); ierr = TSAdaptRegister(TSADAPTCFL, TSAdaptCreate_CFL);CHKERRQ(ierr); ierr = TSAdaptRegister(TSADAPTGLEE, TSAdaptCreate_GLEE);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C TSAdaptFinalizePackage - This function destroys everything in the TS package. It is called from PetscFinalize(). Level: developer .keywords: Petsc, destroy, package .seealso: PetscFinalize() @*/ PetscErrorCode TSAdaptFinalizePackage(void) { PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscFunctionListDestroy(&TSAdaptList);CHKERRQ(ierr); TSAdaptPackageInitialized = PETSC_FALSE; TSAdaptRegisterAllCalled = PETSC_FALSE; PetscFunctionReturn(0); } /*@C TSAdaptInitializePackage - This function initializes everything in the TSAdapt package. It is called from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to TSAdaptCreate() when using static libraries. Level: developer .keywords: TSAdapt, initialize, package .seealso: PetscInitialize() @*/ PetscErrorCode TSAdaptInitializePackage(void) { PetscErrorCode ierr; PetscFunctionBegin; if (TSAdaptPackageInitialized) PetscFunctionReturn(0); TSAdaptPackageInitialized = PETSC_TRUE; ierr = PetscClassIdRegister("TSAdapt",&TSADAPT_CLASSID);CHKERRQ(ierr); ierr = TSAdaptRegisterAll();CHKERRQ(ierr); ierr = PetscRegisterFinalize(TSAdaptFinalizePackage);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C TSAdaptSetType - sets the approach used for the error adapter, currently there is only TSADAPTBASIC and TSADAPTNONE Logicially Collective on TSAdapt Input Parameter: + adapt - the TS adapter, most likely obtained with TSGetAdapt() - type - either TSADAPTBASIC or TSADAPTNONE Options Database: . -ts_adapt_type - to set the adapter type Level: intermediate .keywords: TSAdapt, create .seealso: TSGetAdapt(), TSAdaptDestroy(), TSAdaptType, TSAdaptGetType() @*/ PetscErrorCode TSAdaptSetType(TSAdapt adapt,TSAdaptType type) { PetscBool match; PetscErrorCode ierr,(*r)(TSAdapt); PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidCharPointer(type,2); ierr = PetscObjectTypeCompare((PetscObject)adapt,type,&match);CHKERRQ(ierr); if (match) PetscFunctionReturn(0); ierr = PetscFunctionListFind(TSAdaptList,type,&r);CHKERRQ(ierr); if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown TSAdapt type \"%s\" given",type); if (adapt->ops->destroy) {ierr = (*adapt->ops->destroy)(adapt);CHKERRQ(ierr);} ierr = PetscMemzero(adapt->ops,sizeof(struct _TSAdaptOps));CHKERRQ(ierr); ierr = PetscObjectChangeTypeName((PetscObject)adapt,type);CHKERRQ(ierr); ierr = (*r)(adapt);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C TSAdaptGetType - gets the TS adapter method type (as a string). Not Collective Input Parameter: . adapt - The TS adapter, most likely obtained with TSGetAdapt() Output Parameter: . type - The name of TS adapter method Level: intermediate .keywords: TSAdapt, get, type .seealso TSAdaptSetType() @*/ PetscErrorCode TSAdaptGetType(TSAdapt adapt,TSAdaptType *type) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidPointer(type,2); *type = ((PetscObject)adapt)->type_name; PetscFunctionReturn(0); } PetscErrorCode TSAdaptSetOptionsPrefix(TSAdapt adapt,const char prefix[]) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); ierr = PetscObjectSetOptionsPrefix((PetscObject)adapt,prefix);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C TSAdaptLoad - Loads a TSAdapt that has been stored in binary with TSAdaptView(). Collective on PetscViewer Input Parameters: + newdm - the newly loaded TSAdapt, this needs to have been created with TSAdaptCreate() or some related function before a call to TSAdaptLoad(). - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or HDF5 file viewer, obtained from PetscViewerHDF5Open() Level: intermediate Notes: The type is determined by the data in the file, any type set into the TSAdapt before this call is ignored. Notes for advanced users: Most users should not need to know the details of the binary storage format, since TSAdaptLoad() and TSAdaptView() completely hide these details. But for anyone who's interested, the standard binary matrix storage format is .vb has not yet been determined .ve .seealso: PetscViewerBinaryOpen(), TSAdaptView(), MatLoad(), VecLoad() @*/ PetscErrorCode TSAdaptLoad(TSAdapt adapt,PetscViewer viewer) { PetscErrorCode ierr; PetscBool isbinary; char type[256]; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); ierr = PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);CHKERRQ(ierr); ierr = TSAdaptSetType(adapt,type);CHKERRQ(ierr); if (adapt->ops->load) { ierr = (*adapt->ops->load)(adapt,viewer);CHKERRQ(ierr); } PetscFunctionReturn(0); } PetscErrorCode TSAdaptView(TSAdapt adapt,PetscViewer viewer) { PetscErrorCode ierr; PetscBool iascii,isbinary,isnone; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (!viewer) {ierr = PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)adapt),&viewer);CHKERRQ(ierr);} PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); PetscCheckSameComm(adapt,1,viewer,2); ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); if (iascii) { ierr = PetscObjectPrintClassNamePrefixType((PetscObject)adapt,viewer);CHKERRQ(ierr); ierr = PetscObjectTypeCompare((PetscObject)adapt,TSADAPTNONE,&isnone);CHKERRQ(ierr); if (!isnone) { if (adapt->always_accept) {ierr = PetscViewerASCIIPrintf(viewer," always accepting steps\n");CHKERRQ(ierr);} ierr = PetscViewerASCIIPrintf(viewer," safety factor %g\n",(double)adapt->safety);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," extra safety factor after step rejection %g\n",(double)adapt->reject_safety);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," clip fastest increase %g\n",(double)adapt->clip[1]);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," clip fastest decrease %g\n",(double)adapt->clip[0]);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," maximum allowed timestep %g\n",(double)adapt->dt_max);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer," minimum allowed timestep %g\n",(double)adapt->dt_min);CHKERRQ(ierr); } if (adapt->ops->view) { ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr); ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); } } else if (isbinary) { char type[256]; /* need to save FILE_CLASS_ID for adapt class */ ierr = PetscStrncpy(type,((PetscObject)adapt)->type_name,256);CHKERRQ(ierr); ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); } else if (adapt->ops->view) { ierr = (*adapt->ops->view)(adapt,viewer);CHKERRQ(ierr); } PetscFunctionReturn(0); } /*@ TSAdaptReset - Resets a TSAdapt context. Collective on TS Input Parameter: . adapt - the TSAdapt context obtained from TSAdaptCreate() Level: developer .seealso: TSAdaptCreate(), TSAdaptDestroy() @*/ PetscErrorCode TSAdaptReset(TSAdapt adapt) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (adapt->ops->reset) {ierr = (*adapt->ops->reset)(adapt);CHKERRQ(ierr);} PetscFunctionReturn(0); } PetscErrorCode TSAdaptDestroy(TSAdapt *adapt) { PetscErrorCode ierr; PetscFunctionBegin; if (!*adapt) PetscFunctionReturn(0); PetscValidHeaderSpecific(*adapt,TSADAPT_CLASSID,1); if (--((PetscObject)(*adapt))->refct > 0) {*adapt = NULL; PetscFunctionReturn(0);} ierr = TSAdaptReset(*adapt);CHKERRQ(ierr); if ((*adapt)->ops->destroy) {ierr = (*(*adapt)->ops->destroy)(*adapt);CHKERRQ(ierr);} ierr = PetscViewerDestroy(&(*adapt)->monitor);CHKERRQ(ierr); ierr = PetscHeaderDestroy(adapt);CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ TSAdaptSetMonitor - Monitor the choices made by the adaptive controller Collective on TSAdapt Input Arguments: + adapt - adaptive controller context - flg - PETSC_TRUE to active a monitor, PETSC_FALSE to disable Options Database Keys: . -ts_adapt_monitor - to turn on monitoring Level: intermediate .seealso: TSAdaptChoose() @*/ PetscErrorCode TSAdaptSetMonitor(TSAdapt adapt,PetscBool flg) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidLogicalCollectiveBool(adapt,flg,2); if (flg) { if (!adapt->monitor) {ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)adapt),"stdout",&adapt->monitor);CHKERRQ(ierr);} } else { ierr = PetscViewerDestroy(&adapt->monitor);CHKERRQ(ierr); } PetscFunctionReturn(0); } /*@C TSAdaptSetCheckStage - Set a callback to check convergence for a stage Logically collective on TSAdapt Input Arguments: + adapt - adaptive controller context - func - stage check function Arguments of func: $ PetscErrorCode func(TSAdapt adapt,TS ts,PetscBool *accept) + adapt - adaptive controller context . ts - time stepping context - accept - pending choice of whether to accept, can be modified by this routine Level: advanced .seealso: TSAdaptChoose() @*/ PetscErrorCode TSAdaptSetCheckStage(TSAdapt adapt,PetscErrorCode (*func)(TSAdapt,TS,PetscReal,Vec,PetscBool*)) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); adapt->checkstage = func; PetscFunctionReturn(0); } /*@ TSAdaptSetAlwaysAccept - Set whether to always accept steps regardless of any error or stability condition not meeting the prescribed goal. Logically collective on TSAdapt Input Arguments: + adapt - time step adaptivity context, usually gotten with TSGetAdapt() - flag - whether to always accept steps Options Database Keys: . -ts_adapt_always_accept - to always accept steps Level: intermediate .seealso: TSAdapt, TSAdaptChoose() @*/ PetscErrorCode TSAdaptSetAlwaysAccept(TSAdapt adapt,PetscBool flag) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidLogicalCollectiveBool(adapt,flag,2); adapt->always_accept = flag; PetscFunctionReturn(0); } /*@ TSAdaptSetSafety - Set safety factors Logically collective on TSAdapt Input Arguments: + adapt - adaptive controller context . safety - safety factor relative to target error/stability goal - reject_safety - extra safety factor to apply if the last step was rejected Options Database Keys: + -ts_adapt_safety - to set safety factor - -ts_adapt_reject_safety - to set reject safety factor Level: intermediate .seealso: TSAdapt, TSAdaptGetSafety(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptSetSafety(TSAdapt adapt,PetscReal safety,PetscReal reject_safety) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidLogicalCollectiveReal(adapt,safety,2); PetscValidLogicalCollectiveReal(adapt,reject_safety,3); if (safety != PETSC_DEFAULT && safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be non negative",(double)safety); if (safety != PETSC_DEFAULT && safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Safety factor %g must be less than one",(double)safety); if (reject_safety != PETSC_DEFAULT && reject_safety < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reject safety factor %g must be non negative",(double)reject_safety); if (reject_safety != PETSC_DEFAULT && reject_safety > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reject safety factor %g must be less than one",(double)reject_safety); if (safety != PETSC_DEFAULT) adapt->safety = safety; if (reject_safety != PETSC_DEFAULT) adapt->reject_safety = reject_safety; PetscFunctionReturn(0); } /*@ TSAdaptGetSafety - Get safety factors Not Collective Input Arguments: . adapt - adaptive controller context Ouput Arguments: . safety - safety factor relative to target error/stability goal + reject_safety - extra safety factor to apply if the last step was rejected Level: intermediate .seealso: TSAdapt, TSAdaptSetSafety(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptGetSafety(TSAdapt adapt,PetscReal *safety,PetscReal *reject_safety) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (safety) PetscValidRealPointer(safety,2); if (reject_safety) PetscValidRealPointer(reject_safety,3); if (safety) *safety = adapt->safety; if (reject_safety) *reject_safety = adapt->reject_safety; PetscFunctionReturn(0); } /*@ TSAdaptSetClip - Sets the admissible decrease/increase factor in step size Logically collective on TSAdapt Input Arguments: + adapt - adaptive controller context . low - admissible decrease factor - high - admissible increase factor Options Database Keys: . -ts_adapt_clip , - to set admissible time step decrease and increase factors Level: intermediate .seealso: TSAdaptChoose(), TSAdaptGetClip() @*/ PetscErrorCode TSAdaptSetClip(TSAdapt adapt,PetscReal low,PetscReal high) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidLogicalCollectiveReal(adapt,low,2); PetscValidLogicalCollectiveReal(adapt,high,3); if (low != PETSC_DEFAULT && low < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be non negative",(double)low); if (low != PETSC_DEFAULT && low > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Decrease factor %g must be less than one",(double)low); if (high != PETSC_DEFAULT && high < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Increase factor %g must be geather than one",(double)high); if (low != PETSC_DEFAULT) adapt->clip[0] = low; if (high != PETSC_DEFAULT) adapt->clip[1] = high; PetscFunctionReturn(0); } /*@ TSAdaptGetClip - Gets the admissible decrease/increase factor in step size Not Collective Input Arguments: . adapt - adaptive controller context Ouput Arguments: + low - optional, admissible decrease factor - high - optional, admissible increase factor Level: intermediate .seealso: TSAdaptChoose(), TSAdaptSetClip() @*/ PetscErrorCode TSAdaptGetClip(TSAdapt adapt,PetscReal *low,PetscReal *high) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (low) PetscValidRealPointer(low,2); if (high) PetscValidRealPointer(high,3); if (low) *low = adapt->clip[0]; if (high) *high = adapt->clip[1]; PetscFunctionReturn(0); } /*@ TSAdaptSetStepLimits - Set the minimum and maximum step sizes to be considered by the controller Logically collective on TSAdapt Input Arguments: + adapt - time step adaptivity context, usually gotten with TSGetAdapt() . hmin - minimum time step - hmax - maximum time step Options Database Keys: + -ts_adapt_dt_min - to set minimum time step - -ts_adapt_dt_max - to set maximum time step Level: intermediate .seealso: TSAdapt, TSAdaptGetStepLimits(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptSetStepLimits(TSAdapt adapt,PetscReal hmin,PetscReal hmax) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidLogicalCollectiveReal(adapt,hmin,2); PetscValidLogicalCollectiveReal(adapt,hmax,3); if (hmin != PETSC_DEFAULT && hmin < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmin); if (hmax != PETSC_DEFAULT && hmax < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Minimum time step %g must be non negative",(double)hmax); if (hmin != PETSC_DEFAULT) adapt->dt_min = hmin; if (hmax != PETSC_DEFAULT) adapt->dt_max = hmax; hmin = adapt->dt_min; hmax = adapt->dt_max; if (hmax <= hmin) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Maximum time step %g must geather than minimum time step %g",(double)hmax,(double)hmin); PetscFunctionReturn(0); } /*@ TSAdaptGetStepLimits - Get the minimum and maximum step sizes to be considered by the controller Not Collective Input Arguments: . adapt - time step adaptivity context, usually gotten with TSGetAdapt() Output Arguments: + hmin - minimum time step - hmax - maximum time step Level: intermediate .seealso: TSAdapt, TSAdaptSetStepLimits(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptGetStepLimits(TSAdapt adapt,PetscReal *hmin,PetscReal *hmax) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (hmin) PetscValidRealPointer(hmin,2); if (hmax) PetscValidRealPointer(hmax,3); if (hmin) *hmin = adapt->dt_min; if (hmax) *hmax = adapt->dt_max; PetscFunctionReturn(0); } /* TSAdaptSetFromOptions - Sets various TSAdapt parameters from user options. Collective on TSAdapt Input Parameter: . adapt - the TSAdapt context Options Database Keys: + -ts_adapt_type - algorithm to use for adaptivity . -ts_adapt_always_accept - always accept steps regardless of error/stability goals . -ts_adapt_safety - safety factor relative to target error/stability goal . -ts_adapt_reject_safety - extra safety factor to apply if the last step was rejected . -ts_adapt_clip - admissible time step decrease and increase factors . -ts_adapt_dt_min - minimum timestep to use . -ts_adapt_dt_max - maximum timestep to use . -ts_adapt_scale_solve_failed - scale timestep by this factor if a solve fails - -ts_adapt_wnormtype <2 or infinity> - type of norm for computing error estimates Level: advanced Notes: This function is automatically called by TSSetFromOptions() .keywords: TS, TSGetAdapt(), TSAdaptSetType(), TSAdaptSetStepLimits() .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptSetAlwaysAccept(), TSAdaptSetSafety(), TSAdaptSetClip(), TSAdaptSetStepLimits(), TSAdaptSetMonitor() */ PetscErrorCode TSAdaptSetFromOptions(PetscOptionItems *PetscOptionsObject,TSAdapt adapt) { PetscErrorCode ierr; char type[256] = TSADAPTBASIC; PetscReal safety,reject_safety,clip[2],hmin,hmax; PetscBool set,flg; PetscInt two; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); /* This should use PetscOptionsBegin() if/when this becomes an object used outside of TS, but currently this * function can only be called from inside TSSetFromOptions() */ ierr = PetscOptionsHead(PetscOptionsObject,"TS Adaptivity options");CHKERRQ(ierr); ierr = PetscOptionsFList("-ts_adapt_type","Algorithm to use for adaptivity","TSAdaptSetType",TSAdaptList,((PetscObject)adapt)->type_name ? ((PetscObject)adapt)->type_name : type,type,sizeof(type),&flg);CHKERRQ(ierr); if (flg || !((PetscObject)adapt)->type_name) { ierr = TSAdaptSetType(adapt,type);CHKERRQ(ierr); } ierr = PetscOptionsBool("-ts_adapt_always_accept","Always accept the step","TSAdaptSetAlwaysAccept",adapt->always_accept,&flg,&set);CHKERRQ(ierr); if (set) {ierr = TSAdaptSetAlwaysAccept(adapt,flg);CHKERRQ(ierr);} safety = adapt->safety; reject_safety = adapt->reject_safety; ierr = PetscOptionsReal("-ts_adapt_safety","Safety factor relative to target error/stability goal","TSAdaptSetSafety",safety,&safety,&set);CHKERRQ(ierr); ierr = PetscOptionsReal("-ts_adapt_reject_safety","Extra safety factor to apply if the last step was rejected","TSAdaptSetSafety",reject_safety,&reject_safety,&flg);CHKERRQ(ierr); if (set || flg) {ierr = TSAdaptSetSafety(adapt,safety,reject_safety);CHKERRQ(ierr);} two = 2; clip[0] = adapt->clip[0]; clip[1] = adapt->clip[1]; ierr = PetscOptionsRealArray("-ts_adapt_clip","Admissible decrease/increase factor in step size","TSAdaptSetClip",clip,&two,&set);CHKERRQ(ierr); if (set && (two != 2)) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Must give exactly two values to -ts_adapt_clip"); if (set) {ierr = TSAdaptSetClip(adapt,clip[0],clip[1]);CHKERRQ(ierr);} hmin = adapt->dt_min; hmax = adapt->dt_max; ierr = PetscOptionsReal("-ts_adapt_dt_min","Minimum time step considered","TSAdaptSetStepLimits",hmin,&hmin,&set);CHKERRQ(ierr); ierr = PetscOptionsReal("-ts_adapt_dt_max","Maximum time step considered","TSAdaptSetStepLimits",hmax,&hmax,&flg);CHKERRQ(ierr); if (set || flg) {ierr = TSAdaptSetStepLimits(adapt,hmin,hmax);CHKERRQ(ierr);} ierr = PetscOptionsReal("-ts_adapt_scale_solve_failed","Scale step by this factor if solve fails","",adapt->scale_solve_failed,&adapt->scale_solve_failed,NULL);CHKERRQ(ierr); ierr = PetscOptionsEnum("-ts_adapt_wnormtype","Type of norm computed for error estimation","",NormTypes,(PetscEnum)adapt->wnormtype,(PetscEnum*)&adapt->wnormtype,NULL);CHKERRQ(ierr); if (adapt->wnormtype != NORM_2 && adapt->wnormtype != NORM_INFINITY) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_SUP,"Only 2-norm and infinite norm supported"); ierr = PetscOptionsBool("-ts_adapt_monitor","Print choices made by adaptive controller","TSAdaptSetMonitor",adapt->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr); if (set) {ierr = TSAdaptSetMonitor(adapt,flg);CHKERRQ(ierr);} if (adapt->ops->setfromoptions) {ierr = (*adapt->ops->setfromoptions)(PetscOptionsObject,adapt);CHKERRQ(ierr);} ierr = PetscOptionsTail();CHKERRQ(ierr); PetscFunctionReturn(0); } /*@ TSAdaptCandidatesClear - clear any previously set candidate schemes Logically collective on TSAdapt Input Argument: . adapt - adaptive controller Level: developer .seealso: TSAdapt, TSAdaptCreate(), TSAdaptCandidateAdd(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptCandidatesClear(TSAdapt adapt) { PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); ierr = PetscMemzero(&adapt->candidates,sizeof(adapt->candidates));CHKERRQ(ierr); PetscFunctionReturn(0); } /*@C TSAdaptCandidateAdd - add a candidate scheme for the adaptive controller to select from Logically collective on TSAdapt Input Arguments: + adapt - time step adaptivity context, obtained with TSGetAdapt() or TSAdaptCreate() . name - name of the candidate scheme to add . order - order of the candidate scheme . stageorder - stage order of the candidate scheme . ccfl - stability coefficient relative to explicit Euler, used for CFL constraints . cost - relative measure of the amount of work required for the candidate scheme - inuse - indicates that this scheme is the one currently in use, this flag can only be set for one scheme Note: This routine is not available in Fortran. Level: developer .seealso: TSAdaptCandidatesClear(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptCandidateAdd(TSAdapt adapt,const char name[],PetscInt order,PetscInt stageorder,PetscReal ccfl,PetscReal cost,PetscBool inuse) { PetscInt c; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (order < 1) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Classical order %D must be a positive integer",order); if (inuse) { if (adapt->candidates.inuse_set) SETERRQ(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_WRONGSTATE,"Cannot set the inuse method twice, maybe forgot to call TSAdaptCandidatesClear()"); adapt->candidates.inuse_set = PETSC_TRUE; } /* first slot if this is the current scheme, otherwise the next available slot */ c = inuse ? 0 : !adapt->candidates.inuse_set + adapt->candidates.n; adapt->candidates.name[c] = name; adapt->candidates.order[c] = order; adapt->candidates.stageorder[c] = stageorder; adapt->candidates.ccfl[c] = ccfl; adapt->candidates.cost[c] = cost; adapt->candidates.n++; PetscFunctionReturn(0); } /*@C TSAdaptCandidatesGet - Get the list of candidate orders of accuracy and cost Not Collective Input Arguments: . adapt - time step adaptivity context Output Arguments: + n - number of candidate schemes, always at least 1 . order - the order of each candidate scheme . stageorder - the stage order of each candidate scheme . ccfl - the CFL coefficient of each scheme - cost - the relative cost of each scheme Level: developer Note: The current scheme is always returned in the first slot .seealso: TSAdaptCandidatesClear(), TSAdaptCandidateAdd(), TSAdaptChoose() @*/ PetscErrorCode TSAdaptCandidatesGet(TSAdapt adapt,PetscInt *n,const PetscInt **order,const PetscInt **stageorder,const PetscReal **ccfl,const PetscReal **cost) { PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); if (n) *n = adapt->candidates.n; if (order) *order = adapt->candidates.order; if (stageorder) *stageorder = adapt->candidates.stageorder; if (ccfl) *ccfl = adapt->candidates.ccfl; if (cost) *cost = adapt->candidates.cost; PetscFunctionReturn(0); } /*@C TSAdaptChoose - choose which method and step size to use for the next step Collective on TSAdapt Input Arguments: + adapt - adaptive contoller - h - current step size Output Arguments: + next_sc - optional, scheme to use for the next step . next_h - step size to use for the next step - accept - PETSC_TRUE to accept the current step, PETSC_FALSE to repeat the current step with the new step size Note: The input value of parameter accept is retained from the last time step, so it will be PETSC_FALSE if the step is being retried after an initial rejection. Level: developer .seealso: TSAdapt, TSAdaptCandidatesClear(), TSAdaptCandidateAdd() @*/ PetscErrorCode TSAdaptChoose(TSAdapt adapt,TS ts,PetscReal h,PetscInt *next_sc,PetscReal *next_h,PetscBool *accept) { PetscErrorCode ierr; PetscInt ncandidates = adapt->candidates.n; PetscInt scheme = 0; PetscReal wlte = -1.0; PetscReal wltea = -1.0; PetscReal wlter = -1.0; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidHeaderSpecific(ts,TS_CLASSID,2); if (next_sc) PetscValidIntPointer(next_sc,4); PetscValidPointer(next_h,5); PetscValidIntPointer(accept,6); if (next_sc) *next_sc = 0; /* Do not mess with adaptivity while handling events*/ if (ts->event && ts->event->status != TSEVENT_NONE) { *next_h = h; *accept = PETSC_TRUE; PetscFunctionReturn(0); } ierr = (*adapt->ops->choose)(adapt,ts,h,&scheme,next_h,accept,&wlte,&wltea,&wlter);CHKERRQ(ierr); if (scheme < 0 || (ncandidates > 0 && scheme >= ncandidates)) SETERRQ2(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Chosen scheme %D not in valid range 0..%D",scheme,ncandidates-1); if (*next_h < 0) SETERRQ1(PetscObjectComm((PetscObject)adapt),PETSC_ERR_ARG_OUTOFRANGE,"Computed step size %g must be positive",(double)*next_h); if (next_sc) *next_sc = scheme; if (*accept && ts->exact_final_time == TS_EXACTFINALTIME_MATCHSTEP) { /* Increase/reduce step size if end time of next step is close to or overshoots max time */ PetscReal t = ts->ptime + ts->time_step, h = *next_h; PetscReal tend = t + h, tmax = ts->max_time, hmax = tmax - t; PetscReal a = (PetscReal)1.01; /* allow 1% step size increase */ if (t < tmax && tend > tmax) *next_h = hmax; if (t < tmax && tend < tmax && h > hmax/2) *next_h = hmax/2; if (t < tmax && tend < tmax && h*a > hmax) *next_h = hmax; } if (adapt->monitor) { const char *sc_name = (scheme < ncandidates) ? adapt->candidates.name[scheme] : ""; ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); if (wlte < 0) { ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h);CHKERRQ(ierr); } else { ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s %s %D:%s step %3D %s t=%-11g+%10.3e dt=%-10.3e wlte=%5.3g wltea=%5.3g wlter=%5.3g\n",((PetscObject)adapt)->type_name,((PetscObject)ts)->type_name,scheme,sc_name,ts->steps,*accept ? "accepted" : "rejected",(double)ts->ptime,(double)h,(double)*next_h,(double)wlte,(double)wltea,(double)wlter);CHKERRQ(ierr); } ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); } PetscFunctionReturn(0); } /*@ TSAdaptCheckStage - checks whether to accept a stage, (e.g. reject and change time step size if nonlinear solve fails) Collective on TSAdapt Input Arguments: + adapt - adaptive controller context . ts - time stepper . t - Current simulation time - Y - Current solution vector Output Arguments: . accept - PETSC_TRUE to accept the stage, PETSC_FALSE to reject Level: developer .seealso: @*/ PetscErrorCode TSAdaptCheckStage(TSAdapt adapt,TS ts,PetscReal t,Vec Y,PetscBool *accept) { PetscErrorCode ierr; SNESConvergedReason snesreason = SNES_CONVERGED_ITERATING; PetscFunctionBegin; PetscValidHeaderSpecific(adapt,TSADAPT_CLASSID,1); PetscValidHeaderSpecific(ts,TS_CLASSID,2); PetscValidIntPointer(accept,3); if (ts->snes) {ierr = SNESGetConvergedReason(ts->snes,&snesreason);CHKERRQ(ierr);} if (snesreason < 0) { *accept = PETSC_FALSE; if (++ts->num_snes_failures >= ts->max_snes_failures && ts->max_snes_failures > 0) { ts->reason = TS_DIVERGED_NONLINEAR_SOLVE; ierr = PetscInfo2(ts,"Step=%D, nonlinear solve failures %D greater than current TS allowed, stopping solve\n",ts->steps,ts->num_snes_failures);CHKERRQ(ierr); if (adapt->monitor) { ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected t=%-11g+%10.3e, nonlinear solve failures %D greater than current TS allowed\n",((PetscObject)adapt)->type_name,ts->steps,(double)ts->ptime,(double)ts->time_step,ts->num_snes_failures);CHKERRQ(ierr); ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); } } } else { *accept = PETSC_TRUE; ierr = TSFunctionDomainError(ts,t,Y,accept);CHKERRQ(ierr); if(*accept && adapt->checkstage) { ierr = (*adapt->checkstage)(adapt,ts,t,Y,accept);CHKERRQ(ierr); } } if(!(*accept) && !ts->reason) { PetscReal dt,new_dt; ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr); new_dt = dt * adapt->scale_solve_failed; ierr = TSSetTimeStep(ts,new_dt);CHKERRQ(ierr); adapt->timestepjustincreased += 4; if (adapt->monitor) { ierr = PetscViewerASCIIAddTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(adapt->monitor," TSAdapt %s step %3D stage rejected (%s) t=%-11g+%10.3e retrying with dt=%-10.3e\n",((PetscObject)adapt)->type_name,ts->steps,SNESConvergedReasons[snesreason],(double)ts->ptime,(double)dt,(double)new_dt);CHKERRQ(ierr); ierr = PetscViewerASCIISubtractTab(adapt->monitor,((PetscObject)adapt)->tablevel);CHKERRQ(ierr); } } PetscFunctionReturn(0); } /*@ TSAdaptCreate - create an adaptive controller context for time stepping Collective on MPI_Comm Input Parameter: . comm - The communicator Output Parameter: . adapt - new TSAdapt object Level: developer Notes: TSAdapt creation is handled by TS, so users should not need to call this function. .keywords: TSAdapt, create .seealso: TSGetAdapt(), TSAdaptSetType(), TSAdaptDestroy() @*/ PetscErrorCode TSAdaptCreate(MPI_Comm comm,TSAdapt *inadapt) { PetscErrorCode ierr; TSAdapt adapt; PetscFunctionBegin; PetscValidPointer(inadapt,1); *inadapt = NULL; ierr = TSAdaptInitializePackage();CHKERRQ(ierr); ierr = PetscHeaderCreate(adapt,TSADAPT_CLASSID,"TSAdapt","Time stepping adaptivity","TS",comm,TSAdaptDestroy,TSAdaptView);CHKERRQ(ierr); adapt->always_accept = PETSC_FALSE; adapt->safety = 0.9; adapt->reject_safety = 0.5; adapt->clip[0] = 0.1; adapt->clip[1] = 10.; adapt->dt_min = 1e-20; adapt->dt_max = 1e+20; adapt->scale_solve_failed = 0.25; adapt->wnormtype = NORM_2; *inadapt = adapt; PetscFunctionReturn(0); }