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 namespace Tempus {
20 
57 template <class Scalar>
59  : virtual public TimeStepControlStrategy<Scalar> {
60  public:
63  : controller_("PID"),
64  KI_(0.58),
65  KP_(0.21),
66  KD_(0.10),
67  safetyFactor_(0.90),
69  facMax_(5.0),
70  facMin_(0.5)
71  {
73 
74  this->setStrategyType("Integral Controller");
75  this->setStepType("Variable");
76  this->setName("Integral Controller");
77  this->initialize();
78  }
79 
82  std::string controller, Scalar KI, Scalar KP, Scalar KD,
83  Scalar safetyFactor, Scalar safetyFactorAfterReject, Scalar facMax,
84  Scalar facMin, std::string name = "Integral Controller")
85  : controller_(controller),
86  KI_(KI),
87  KP_(KP),
88  KD_(KD),
89  safetyFactor_(safetyFactor),
90  safetyFactorAfterReject_(safetyFactorAfterReject),
91  facMax_(facMax),
92  facMin_(facMin)
93  {
95 
96  this->setStrategyType("Integral Controller");
97  this->setStepType("Variable");
98  this->setName(name);
99  this->initialize();
100  }
101 
104 
106  virtual void setNextTimeStep(
107  const TimeStepControl<Scalar> &tsc,
108  Teuchos::RCP<SolutionHistory<Scalar> > solutionHistory,
109  Status & /* integratorStatus */) override
110  {
111  using Teuchos::RCP;
112 
113  this->checkInitialized();
114 
115  // Take first step with initial time step provided
116  if (!firstSuccessfulStep_) {
117  firstSuccessfulStep_ = true;
118  return;
119  }
120 
121  RCP<SolutionState<Scalar> > workingState =
122  solutionHistory->getWorkingState();
123  Scalar beta = 1.0;
124 
125  // assumes the embedded solution is the low order solution
126  int order = workingState->getOrder() - 1;
127  Scalar dt = workingState->getTimeStep();
128 
129  // Get the relative errors.
130  Scalar errN = workingState->getErrorRel();
131  Scalar errNm1 = workingState->getErrorRelNm1();
132  Scalar errNm2 = workingState->getErrorRelNm2();
133 
134  if (errN < numericalTol<Scalar>()) errN = 1.0;
135  if (errNm1 < numericalTol<Scalar>()) errNm1 = 1.0;
136  if (errNm2 < numericalTol<Scalar>()) errNm2 = 1.0;
137 
138  Scalar k1 = Teuchos::as<Scalar>(-KI_ / order);
139  Scalar k2 = Teuchos::as<Scalar>(KP_ / order);
140  Scalar k3 = Teuchos::as<Scalar>(-KD_ / order);
141 
142  k1 = std::pow(errN, k1);
143  k2 = std::pow(errNm1, k2);
144  k3 = std::pow(errNm2, k3);
145 
146  if (controller_ == "I")
147  beta = safetyFactor_ * k1;
148  else if (controller_ == "PI")
149  beta = safetyFactor_ * k1 * k2;
150  else // (controller_ == "PID")
151  beta = safetyFactor_ * k1 * k2 * k3;
152 
153  beta = std::max(facMin_, beta);
154  beta = std::min(facMax_, beta);
155 
156  // new (optimal) suggested time step
157  dt = beta * dt;
158 
159  if (workingState->getSolutionStatus() == Status::PASSED ||
160  workingState->getSolutionStatus() == Status::WORKING) {
161  if (lastStepRejected_) {
162  dt = std::min(dt, workingState->getTimeStep());
163  }
164  else {
166  }
167  lastStepRejected_ = false;
168  }
169  else {
171  lastStepRejected_ = true;
172  }
173 
174  // update dt
175  workingState->setTimeStep(dt);
176  workingState->setTime(solutionHistory->getCurrentState()->getTime() + dt);
177  }
178 
180 
181  std::string description() const override
182  {
183  return "Tempus::TimeStepControlStrategyIntegralController";
184  }
185 
187  const Teuchos::EVerbosityLevel verbLevel) const override
188  {
189  auto l_out = Teuchos::fancyOStream(out.getOStream());
190  Teuchos::OSTab ostab(*l_out, 2, this->description());
191  l_out->setOutputToRootOnly(0);
192 
193  *l_out << "\n--- " << this->description() << " ---" << std::endl;
194 
195  if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_MEDIUM)) {
196  *l_out << " Strategy Type = "
197  << this->getStrategyType() << std::endl
198  << " Step Type = " << this->getStepType()
199  << std::endl
200  << " Controller Type = " << getController()
201  << std::endl
202  << " KI = " << getKI()
203  << std::endl
204  << " KP = " << getKP()
205  << std::endl
206  << " KD = " << getKD()
207  << std::endl
208  << " Safety Factor = " << getSafetyFactor()
209  << std::endl
210  << " Safety Factor After Step Rejection = "
211  << getSafetyFactorAfterReject() << std::endl
212  << " Maximum Safety Factor (INPUT) = " << facMaxINPUT_
213  << std::endl
214  << " Maximum Safety Factor = " << getFacMax()
215  << std::endl
216  << " Minimum Safety Factor = " << getFacMin()
217  << std::endl;
218  *l_out << std::string(this->description().length() + 8, '-') << std::endl;
219  }
220  }
222 
225  const override
226  {
228  Teuchos::parameterList("Time Step Control Strategy");
229 
230  pl->set<std::string>("Strategy Type", this->getStrategyType(),
231  "Integral Controller");
232  pl->set<std::string>("Controller Type", getController(),
233  "Proportional-Integral-Derivative");
234  pl->set<Scalar>("KI", getKI(), "Integral gain");
235  pl->set<Scalar>("KP", getKP(), "Proportional gain");
236  pl->set<Scalar>("KD", getKD(), "Derivative gain");
237  pl->set<Scalar>("Safety Factor", getSafetyFactor(), "Safety Factor");
238  pl->set<Scalar>("Safety Factor After Step Rejection",
240  "Safety Factor Following Step Rejection");
241  pl->set<Scalar>("Maximum Safety Factor", getFacMax(),
242  "Maximum Safety Factor");
243  pl->set<Scalar>("Minimum Safety Factor", getFacMin(),
244  "Minimum Safety Factor");
245  return pl;
246  }
247 
248  virtual void initialize() const override
249  {
250  TEUCHOS_TEST_FOR_EXCEPTION(safetyFactor_ <= 0.0, std::out_of_range,
251  "Error - Invalid value of Safety Factory= "
252  << safetyFactor_ << "! \n"
253  << "Safety Factor must be > 0.0.\n");
254 
256  facMax_ <= 0.0, std::out_of_range,
257  "Error - Invalid value of Maximum Safety Factory= "
258  << facMax_ << "! \n"
259  << "Maximum Safety Factor must be > 0.0.\n");
260 
262  facMax_ <= 0.0, std::out_of_range,
263  "Error - Invalid value of Minimum Safety Factory= "
264  << facMin_ << "! \n"
265  << "Minimum Safety Factor must be > 0.0.\n");
266 
267  TEUCHOS_TEST_FOR_EXCEPTION(((controller_ != "I") && (controller_ != "PI") &&
268  (controller_ != "PID")),
269  std::invalid_argument,
270  "Error - Invalid choice of Controller Type = "
271  << controller_ << "! \n"
272  << "Valid Choice are ['I', 'PI', 'PID'].\n");
273 
274  this->isInitialized_ = true; // Only place where this is set to true!
275  }
276 
277  virtual std::string getController() const { return controller_; }
278  virtual Scalar getKI() const { return KI_; }
279  virtual Scalar getKP() const { return KP_; }
280  virtual Scalar getKD() const { return KD_; }
281  virtual Scalar getSafetyFactor() const { return safetyFactor_; }
282  virtual Scalar getSafetyFactorAfterReject() const
283  {
285  }
286  virtual Scalar getFacMax() const { return facMax_; }
287  virtual Scalar getFacMin() const { return facMin_; }
288 
289  virtual void setController(std::string c)
290  {
291  controller_ = c;
292  this->isInitialized_ = false;
293  }
294  virtual void setKI(Scalar k)
295  {
296  KI_ = k;
297  this->isInitialized_ = false;
298  }
299  virtual void setKP(Scalar k)
300  {
301  KP_ = k;
302  this->isInitialized_ = false;
303  }
304  virtual void setKD(Scalar k)
305  {
306  KD_ = k;
307  this->isInitialized_ = false;
308  }
309  virtual void setSafetyFactor(Scalar f)
310  {
311  safetyFactor_ = f;
312  this->isInitialized_ = false;
313  }
314  virtual void setSafetyFactorAfterReject(Scalar f)
315  {
317  this->isInitialized_ = false;
318  }
319  virtual void setFacMax(Scalar f)
320  {
321  facMax_ = f;
322  facMaxINPUT_ = f;
323  this->isInitialized_ = false;
324  }
325  virtual void setFacMin(Scalar f)
326  {
327  facMin_ = f;
328  this->isInitialized_ = false;
329  }
330 
331  private:
332  std::string controller_;
333  Scalar KI_;
334  Scalar KP_;
335  Scalar KD_;
336  Scalar safetyFactor_;
338  Scalar facMaxINPUT_;
339  Scalar facMax_;
340  Scalar facMin_;
341  bool firstSuccessfulStep_ = false;
342  bool lastStepRejected_ = false;
343 };
344 
345 // Nonmember constructor.
346 template <class Scalar>
350  std::string name = "Integral Controller")
351 {
352  using Teuchos::rcp;
354  if (pList == Teuchos::null || pList->numParams() == 0) return tscs;
355 
357  pList->get<std::string>("Strategy Type") != "Integral Controller",
358  std::logic_error,
359  "Error - Strategy Type != 'Integral Controller'. (='" +
360  pList->get<std::string>("Strategy Type") + "')\n");
361 
362  pList->validateParametersAndSetDefaults(*tscs->getValidParameters());
363 
364  tscs->setController(pList->get<std::string>("Controller Type"));
365  tscs->setKI(pList->get<Scalar>("KI"));
366  tscs->setKP(pList->get<Scalar>("KP"));
367  tscs->setKD(pList->get<Scalar>("KD"));
368  tscs->setSafetyFactor(pList->get<Scalar>("Safety Factor"));
369  tscs->setSafetyFactorAfterReject(
370  pList->get<Scalar>("Safety Factor After Step Rejection"));
371  tscs->setFacMax(pList->get<Scalar>("Maximum Safety Factor"));
372  tscs->setFacMin(pList->get<Scalar>("Minimum Safety Factor"));
373 
374  tscs->setName(name);
375  tscs->initialize();
376 
377  return tscs;
378 }
379 
381 template <class Scalar>
384 {
386  return Teuchos::rcp_const_cast<Teuchos::ParameterList>(
387  t->getValidParameters());
388 }
389 
390 } // namespace Tempus
391 #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.