Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_TimeStepControlStrategyIntegralController.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ****************************************************************************
3 // Tempus: Copyright (2017) Sandia Corporation
4 //
5 // Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
6 // ****************************************************************************
7 // @HEADER
8 
9 #ifndef Tempus_TimeStepControlStrategy_IntegralController_hpp
10 #define Tempus_TimeStepControlStrategy_IntegralController_hpp
11 
12 #include "Tempus_config.hpp"
15 #include "Tempus_SolutionState.hpp"
16 #include "Tempus_SolutionHistory.hpp"
17 #include "Tempus_StepperState.hpp"
18 
19 
20 namespace Tempus {
21 
58 template<class Scalar>
60  : virtual public TimeStepControlStrategy<Scalar>
61 {
62 public:
63 
66  : controller_("PID"), KI_(0.58), KP_(0.21), KD_(0.10),
68  facMax_(5.0), facMin_(0.5)
69  {
71 
72  this->setStrategyType("Integral Controller");
73  this->setStepType("Variable");
74  this->setName("Integral Controller");
75  this->initialize();
76  }
77 
80  Scalar KI, Scalar KP, Scalar KD,
81  Scalar safetyFactor, Scalar safetyFactorAfterReject,
82  Scalar facMax, Scalar facMin, std::string name = "Integral Controller")
83  : controller_(controller), KI_(KI), KP_(KP), KD_(KD),
84  safetyFactor_(safetyFactor),
85  safetyFactorAfterReject_(safetyFactorAfterReject),
86  facMax_(facMax), facMin_(facMin)
87  {
89 
90  this->setStrategyType("Integral Controller");
91  this->setStepType("Variable");
92  this->setName(name);
93  this->initialize();
94  }
95 
96 
99 
101  virtual void setNextTimeStep(const TimeStepControl<Scalar> & tsc,
102  Teuchos::RCP<SolutionHistory<Scalar> > solutionHistory,
103  Status & /* integratorStatus */) override
104  {
105  using Teuchos::RCP;
106 
107  this->checkInitialized();
108 
109  // Take first step with initial time step provided
110  if (!firstSuccessfulStep_){
111  firstSuccessfulStep_ = true;
112  return;
113  }
114 
115  RCP<SolutionState<Scalar> > workingState=solutionHistory->getWorkingState();
116  Scalar beta = 1.0;
117 
118  // assumes the embedded solution is the low order solution
119  int order = workingState->getOrder() - 1;
120  Scalar dt = workingState->getTimeStep();
121 
122  // Get the relative errors.
123  Scalar errN = workingState->getErrorRel();
124  Scalar errNm1 = workingState->getErrorRelNm1();
125  Scalar errNm2 = workingState->getErrorRelNm2();
126 
127  if ( errN < numericalTol<Scalar>()) errN = 1.0;
128  if ( errNm1 < numericalTol<Scalar>()) errNm1 = 1.0;
129  if ( errNm2 < numericalTol<Scalar>()) errNm2 = 1.0;
130 
131  Scalar k1 = Teuchos::as<Scalar>(-KI_ / order);
132  Scalar k2 = Teuchos::as<Scalar>( KP_ / order);
133  Scalar k3 = Teuchos::as<Scalar>(-KD_ / order);
134 
135  k1 = std::pow(errN, k1);
136  k2 = std::pow(errNm1, k2);
137  k3 = std::pow(errNm2, k3);
138 
139  if (controller_ == "I")
140  beta = safetyFactor_*k1;
141  else if (controller_ == "PI")
142  beta = safetyFactor_ *k1*k2;
143  else // (controller_ == "PID")
144  beta = safetyFactor_*k1*k2*k3;
145 
146  beta = std::max(facMin_, beta);
147  beta = std::min(facMax_, beta);
148 
149  // new (optimal) suggested time step
150  dt = beta * dt;
151 
152  if (workingState->getSolutionStatus() == Status::PASSED ||
153  workingState->getSolutionStatus() == Status::WORKING) {
154  if(lastStepRejected_){
155  dt = std::min(dt, workingState->getTimeStep());
156  } else {
158  }
159  lastStepRejected_ = false;
160  } else {
162  lastStepRejected_ = true;
163  }
164 
165  // update dt
166  workingState->setTimeStep(dt);
167  workingState->setTime(solutionHistory->getCurrentState()->getTime() + dt);
168  }
169 
170 
172 
173  std::string description() const override
174  { return "Tempus::TimeStepControlStrategyIntegralController"; }
175 
177  const Teuchos::EVerbosityLevel verbLevel) const override
178  {
179  auto l_out = Teuchos::fancyOStream( out.getOStream() );
180  Teuchos::OSTab ostab(*l_out, 2, this->description());
181  l_out->setOutputToRootOnly(0);
182 
183  *l_out << "\n--- " << this->description() << " ---" << std::endl;
184 
185  if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
186  *l_out << " Strategy Type = " << this->getStrategyType() << std::endl
187  << " Step Type = " << this->getStepType() << std::endl
188  << " Controller Type = " << getController() << std::endl
189  << " KI = " << getKI() << std::endl
190  << " KP = " << getKP() << std::endl
191  << " KD = " << getKD() << std::endl
192  << " Safety Factor = " << getSafetyFactor() << std::endl
193  << " Safety Factor After Step Rejection = " << getSafetyFactorAfterReject() << std::endl
194  << " Maximum Safety Factor (INPUT) = " << facMaxINPUT_ << std::endl
195  << " Maximum Safety Factor = " << getFacMax() << std::endl
196  << " Minimum Safety Factor = " << getFacMin() << std::endl;
197  *l_out << std::string(this->description().length()+8, '-') <<std::endl;
198  }
199  }
201 
204  {
206  Teuchos::parameterList("Time Step Control Strategy");
207 
208  pl->set<std::string>("Strategy Type", this->getStrategyType(), "Integral Controller");
209  pl->set<std::string>("Controller Type", getController(),
210  "Proportional-Integral-Derivative");
211  pl->set<Scalar>("KI" , getKI(), "Integral gain");
212  pl->set<Scalar>("KP" , getKP(), "Proportional gain");
213  pl->set<Scalar>("KD" , getKD(), "Derivative gain");
214  pl->set<Scalar>("Safety Factor" , getSafetyFactor(), "Safety Factor");
215  pl->set<Scalar>("Safety Factor After Step Rejection",
217  "Safety Factor Following Step Rejection");
218  pl->set<Scalar>("Maximum Safety Factor" , getFacMax(), "Maximum Safety Factor");
219  pl->set<Scalar>("Minimum Safety Factor" , getFacMin(), "Minimum Safety Factor");
220  return pl;
221  }
222 
223 
224  virtual void initialize() const override
225  {
226  TEUCHOS_TEST_FOR_EXCEPTION(safetyFactor_ <= 0.0, std::out_of_range,
227  "Error - Invalid value of Safety Factory= " << safetyFactor_ << "! \n"
228  << "Safety Factor must be > 0.0.\n");
229 
230  TEUCHOS_TEST_FOR_EXCEPTION(facMax_ <= 0.0, std::out_of_range,
231  "Error - Invalid value of Maximum Safety Factory= " << facMax_ << "! \n"
232  << "Maximum Safety Factor must be > 0.0.\n");
233 
234  TEUCHOS_TEST_FOR_EXCEPTION(facMax_<= 0.0, std::out_of_range,
235  "Error - Invalid value of Minimum Safety Factory= " << facMin_ << "! \n"
236  << "Minimum Safety Factor must be > 0.0.\n");
237 
239  (controller_ != "PI") &&
240  (controller_ != "PID")), std::invalid_argument,
241  "Error - Invalid choice of Controller Type = " << controller_ << "! \n"
242  << "Valid Choice are ['I', 'PI', 'PID'].\n");
243 
244  this->isInitialized_ = true; // Only place where this is set to true!
245  }
246 
247 
248  virtual std::string getController() const { return controller_; }
249  virtual Scalar getKI() const { return KI_; }
250  virtual Scalar getKP() const { return KP_; }
251  virtual Scalar getKD() const { return KD_; }
252  virtual Scalar getSafetyFactor() const { return safetyFactor_; }
253  virtual Scalar getSafetyFactorAfterReject() const { return safetyFactorAfterReject_; }
254  virtual Scalar getFacMax() const { return facMax_; }
255  virtual Scalar getFacMin() const { return facMin_; }
256 
257  virtual void setController(std::string c) { controller_ = c; this->isInitialized_ = false; }
258  virtual void setKI(Scalar k) { KI_ = k; this->isInitialized_ = false; }
259  virtual void setKP(Scalar k) { KP_ = k; this->isInitialized_ = false; }
260  virtual void setKD(Scalar k) { KD_ = k; this->isInitialized_ = false; }
261  virtual void setSafetyFactor(Scalar f) { safetyFactor_ = f; this->isInitialized_ = false; }
262  virtual void setSafetyFactorAfterReject(Scalar f) { safetyFactorAfterReject_ = f; this->isInitialized_ = false; }
263  virtual void setFacMax(Scalar f) { facMax_ = f; facMaxINPUT_ = f; this->isInitialized_ = false; }
264  virtual void setFacMin(Scalar f) { facMin_ = f; this->isInitialized_ = false; }
265 
266 private:
267 
268  std::string controller_;
269  Scalar KI_;
270  Scalar KP_;
271  Scalar KD_;
272  Scalar safetyFactor_;
274  Scalar facMaxINPUT_;
275  Scalar facMax_;
276  Scalar facMin_;
277  bool firstSuccessfulStep_ = false;
278  bool lastStepRejected_ = false;
279 
280 };
281 
282 
283 // Nonmember constructor.
284 template <class Scalar>
288  std::string name = "Integral Controller")
289 {
290  using Teuchos::rcp;
292  if (pList == Teuchos::null || pList->numParams() == 0) return tscs;
293 
295  pList->get<std::string>("Strategy Type") !=
296  "Integral Controller", std::logic_error,
297  "Error - Strategy Type != 'Integral Controller'. (='"
298  +pList->get<std::string>("Strategy Type")+"')\n");
299 
300  pList->validateParametersAndSetDefaults(*tscs->getValidParameters());
301 
302  tscs->setController (pList->get<std::string>("Controller Type"));
303  tscs->setKI (pList->get<Scalar>("KI"));
304  tscs->setKP (pList->get<Scalar>("KP"));
305  tscs->setKD (pList->get<Scalar>("KD"));
306  tscs->setSafetyFactor(pList->get<Scalar>("Safety Factor"));
307  tscs->setSafetyFactorAfterReject(pList->get<Scalar>("Safety Factor After Step Rejection"));
308  tscs->setFacMax (pList->get<Scalar>("Maximum Safety Factor"));
309  tscs->setFacMin (pList->get<Scalar>("Minimum Safety Factor"));
310 
311  tscs->setName(name);
312  tscs->initialize();
313 
314  return tscs;
315 }
316 
317 
319 template<class Scalar>
321 {
323  return Teuchos::rcp_const_cast<Teuchos::ParameterList> (t->getValidParameters());
324 }
325 
326 
327 } // namespace Tempus
328 #endif // Tempus_TimeStepControlStrategy_IntegralController_hpp
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Return ParameterList with current values.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Ordinal numParams() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Status
Status for the Integrator, the Stepper and the SolutionState.
bool isInitialized_
Bool if strategy is initialized.
TimeStepControl manages the time step size. There several mechanisms that effect the time step size a...
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const override
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
std::string controller_
Control type [&#39;I&#39;, &#39;PI&#39;, &#39;PID&#39;].
Teuchos::RCP< Teuchos::ParameterList > getTimeStepControlStrategyIntegralControllerPL()
Nonmember function to return ParameterList with default values.
TimeStepControlStrategyIntegralController(std::string controller, Scalar KI, Scalar KP, Scalar KD, Scalar safetyFactor, Scalar safetyFactorAfterReject, Scalar facMax, Scalar facMin, std::string name="Integral Controller")
Full Constructor.
virtual void setNextTimeStep(const TimeStepControl< Scalar > &tsc, Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory, Status &) override
Set the time step size.
Teuchos::RCP< TimeStepControlStrategyIntegralController< Scalar > > createTimeStepControlStrategyIntegralController(const Teuchos::RCP< Teuchos::ParameterList > pList, std::string name="Integral Controller")
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
TimeStepControlStrategy class for TimeStepControl.