Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_TimeEventRange_impl.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_TimeEventRange_impl_hpp
10 #define Tempus_TimeEventRange_impl_hpp
11 
13 
14 namespace Tempus {
15 
16 template <class Scalar>
18  : start_(0.0),
19  stop_(0.0),
20  stride_(0.0),
21  numEvents_(1),
22  timeScale_(1.0),
23  relTol_(this->getDefaultTol()),
24  absTol_(this->getDefaultTol()),
25  landOnExactly_(true)
26 {
27  this->setType("Range");
28  setRelTol(this->getDefaultTol()), setTimeRange(0.0, 0.0, 0.0);
29  setLandOnExactly(true);
30  std::ostringstream oss;
31  oss << "TimeEventRange (" << start_ << "; " << stop_ << "; " << stride_
32  << ")";
33  this->setName(oss.str());
34 }
35 
36 template <class Scalar>
37 TimeEventRange<Scalar>::TimeEventRange(Scalar start, Scalar stop, Scalar stride,
38  std::string name, bool landOnExactly,
39  Scalar relTol)
40  : start_(start),
41  stop_(stop),
42  stride_(stride),
43  numEvents_(0),
44  timeScale_(std::max(std::abs(start_), std::abs(stop_))),
45  relTol_(relTol),
46  absTol_(relTol_ * timeScale_),
47  landOnExactly_(landOnExactly)
48 {
49  this->setType("Range");
50  if (name == "") {
51  std::ostringstream oss;
52  oss << "TimeEventRange (" << start << "; " << stop << "; " << stride << ")";
53  this->setName(oss.str());
54  }
55  else {
56  this->setName(name);
57  }
58  setRelTol(relTol);
59  setTimeRange(start, stop, stride);
60  setLandOnExactly(landOnExactly);
61 }
62 
63 template <class Scalar>
64 TimeEventRange<Scalar>::TimeEventRange(Scalar start, Scalar stop, int numEvents,
65  std::string name, bool landOnExactly,
66  Scalar relTol)
67  : start_(start),
68  stop_(stop),
69  stride_(0.0),
70  numEvents_(numEvents),
71  timeScale_(std::max(std::abs(start_), std::abs(stop_))),
72  relTol_(relTol),
73  absTol_(relTol_ * timeScale_),
74  landOnExactly_(landOnExactly)
75 {
76  if (name == "") {
77  std::ostringstream oss;
78  oss << "TimeEventRange (" << start << "; " << stop << "; " << numEvents
79  << ")";
80  this->setName(oss.str());
81  }
82  else {
83  this->setName(name);
84  }
85  setRelTol(relTol);
86  setTimeRange(start, stop, numEvents);
87  setLandOnExactly(landOnExactly);
88 }
89 
90 template <class Scalar>
91 void TimeEventRange<Scalar>::setTimeRange(Scalar start, Scalar stop,
92  Scalar stride)
93 {
94  start_ = start;
95  stop_ = stop;
96  if (stop_ < start_) {
97  Scalar tmp = start_;
98  start_ = stop_;
99  stop_ = tmp;
100  }
101  setTimeScale();
102  setTimeStride(stride);
103 }
104 
105 template <class Scalar>
106 void TimeEventRange<Scalar>::setTimeRange(Scalar start, Scalar stop,
107  int numEvents)
108 {
109  start_ = start;
110  stop_ = stop;
111  if (stop_ < start_) {
112  Scalar tmp = start_;
113  start_ = stop_;
114  stop_ = tmp;
115  }
116  setTimeScale();
117  setNumEvents(numEvents);
118 }
119 
120 template <class Scalar>
122 {
123  start_ = start;
124  if (stop_ < start_) stop_ = start_;
125  setTimeScale();
126  setTimeStride(stride_); // Reset numEvents with the current stride.
127 }
128 
129 template <class Scalar>
131 {
132  stop_ = stop;
133  if (start_ > stop_) start_ = stop_;
134  setTimeScale();
135  setTimeStride(stride_); // Reset numEvents with the current stride.
136 }
137 
138 template <class Scalar>
140 {
141  timeScale_ = std::max(std::abs(start_), std::abs(stop_));
142  absTol_ = relTol_ * timeScale_;
143 
144  // Check if timeScale is near zero.
145  if (approxZero(timeScale_, absTol_)) {
146  timeScale_ = 1.0;
147  absTol_ = relTol_ * timeScale_;
148  }
149 }
150 
151 template <class Scalar>
153 {
155  if (approxEqualAbsTol(start_, stop_, absTol_)) {
156  stride_ = 0.0;
157  numEvents_ = 1;
158  return;
159  }
160 
161  if ((stride_ > stop_ - start_) || (stride_ < 2 * absTol_)) {
162  stride_ = stop_ - start_;
163  }
164 
165  numEvents_ = int((stop_ + absTol_ - start_) / stride_) + 1;
166 }
167 
168 template <class Scalar>
170 {
171  numEvents_ = numEvents;
172  if (numEvents_ < 0) { // Do not use numEvents_ to set stride! Reset
173  // numEvents_ from stride_.
174  if (stride_ < 2 * absTol_) stride_ = 2 * absTol_;
175  numEvents_ = int((stop_ + absTol_ - start_) / stride_) + 1;
176  return;
177  }
178  else if (approxEqualAbsTol(start_, stop_, absTol_)) {
179  numEvents_ = 1;
180  stride_ = 0.0;
181  stride_ = stop_ - start_;
182  }
183  else {
184  if (numEvents_ < 2) numEvents_ = 2;
185  stride_ = (stop_ - start_) / Scalar(numEvents_ - 1);
186  }
187 
188  // If stride_ is smaller than twice the absolute tolerance,
189  // the time steps cannot be differentiated!
190  if (stride_ <= 2 * absTol_) setTimeStride(2 * absTol_);
191 }
192 
193 template <class Scalar>
195 {
196  relTol_ = std::abs(relTol);
197  setTimeScale();
198 }
199 
200 template <class Scalar>
201 bool TimeEventRange<Scalar>::isTime(Scalar time) const
202 {
203  // Check if before first event.
204  if (time < start_ - absTol_) return false;
205 
206  // Check if after last event.
207  const Scalar timeOfLast = start_ + (numEvents_ - 1) * stride_;
208  if (time > timeOfLast + absTol_) return false;
209 
210  int numStrides = 0;
211  if (!approxZero(stride_, 2 * absTol_)) numStrides = (time - start_) / stride_;
212 
213  numStrides = std::min(std::max(0, numStrides), int(numEvents_ - 1));
214  const Scalar leftBracket = start_ + numStrides * stride_;
215 
216  // Check if close to left bracket.
217  if (approxEqualAbsTol(time, leftBracket, absTol_)) return true;
218 
219  // Check if close to right bracket.
220  const Scalar rightBracket = leftBracket + stride_;
221  if (approxEqualAbsTol(time, rightBracket, absTol_)) return true;
222 
223  return false;
224 }
225 
226 template <class Scalar>
228 {
229  return timeOfNextEvent(time) - time; // Neg. time indicates in the past.
230 }
231 
232 template <class Scalar>
234 {
235  // Check if before first event.
236  if (time < start_ - absTol_) return start_;
237 
238  const Scalar timeOfLast = start_ + (numEvents_ - 1) * stride_;
239  // Check if after or close to last event.
240  if (time > timeOfLast - absTol_) return std::numeric_limits<Scalar>::max();
241 
242  int numStrides = 0;
243  if (!approxZero(stride_, 2 * absTol_))
244  numStrides = int((time - start_) / stride_) + 1;
245  const Scalar timeEvent = start_ + numStrides * stride_;
246 
247  // Check timeEvent is near time. If so, return the next event.
248  if (approxEqualAbsTol(time, timeEvent, absTol_)) return timeEvent + stride_;
249 
250  return timeEvent;
251 }
252 
253 template <class Scalar>
254 bool TimeEventRange<Scalar>::eventInRange(Scalar time1, Scalar time2) const
255 {
256  if (time1 > time2) {
257  Scalar tmp = time1;
258  time1 = time2;
259  time2 = tmp;
260  }
261 
262  // Check if range is completely outside time events.
263  const Scalar timeOfLast = start_ + (numEvents_ - 1) * stride_;
264  if (time2 < start_ - absTol_ || timeOfLast + absTol_ < time1) return false;
265 
266  if (approxZero(stride_))
267  return (time1 < start_ - absTol_ && start_ - absTol_ <= time2);
268 
269  const int strideJustBeforeTime1 = std::min(
270  int(numEvents_ - 1),
271  std::max(int(0), int((time1 - start_ + absTol_) / stride_ - 0.5)));
272 
273  const int strideJustAfterTime2 = std::max(
274  int(0), std::min(int(numEvents_ - 1),
275  int((time2 - start_ + absTol_) / stride_ + 0.5)));
276 
277  for (int i = strideJustBeforeTime1; i <= strideJustAfterTime2; i++) {
278  const Scalar timeEvent = start_ + i * stride_;
279  if (time1 < timeEvent - absTol_ && timeEvent - absTol_ <= time2)
280  return true;
281  }
282 
283  return false;
284 }
285 
286 template <class Scalar>
288  Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
289 {
290  auto l_out = Teuchos::fancyOStream(out.getOStream());
291  Teuchos::OSTab ostab(*l_out, 2, "TimeEventRange");
292  l_out->setOutputToRootOnly(0);
293 
294  *l_out << "TimeEventRange:"
295  << "\n"
296  << " name = " << this->getName() << "\n"
297  << " Type = " << this->getType() << "\n"
298  << " start_ = " << start_ << "\n"
299  << " stop_ = " << stop_ << "\n"
300  << " stride_ = " << stride_ << "\n"
301  << " numEvents_ = " << numEvents_ << "\n"
302  << " timeScale_ = " << timeScale_ << "\n"
303  << " relTol_ = " << relTol_ << "\n"
304  << " absTol_ = " << absTol_ << "\n"
305  << " landOnExactly_ = " << landOnExactly_ << std::endl;
306 }
307 
308 template <class Scalar>
311 {
313  Teuchos::parameterList("Time Event Range");
314 
315  pl->setName(this->getName());
316  pl->set("Name", this->getName());
317  pl->set("Type", this->getType());
318 
319  pl->set("Start Time", getTimeStart(), "Start of time range");
320  pl->set("Stop Time", getTimeStop(), "Stop of time range");
321  pl->set("Stride Time", getTimeStride(), "Stride of time range");
322 
323  if (getTimeStride() * Scalar(getNumEvents() - 1) -
324  (getTimeStop() - getTimeStart()) <
325  getAbsTol())
326  pl->set("Number of Events", getNumEvents(),
327  "Number of events in time range. If specified, 'Stride Time' is "
328  "reset.");
329 
330  pl->set("Relative Tolerance", getRelTol(),
331  "Relative time tolerance for matching time events.");
332 
333  pl->set("Land On Exactly", getLandOnExactly(),
334  "Should these time events be landed on exactly, i.e, adjust the "
335  "timestep to hit time event, versus stepping over and keeping the "
336  "time step unchanged.");
337 
338  return pl;
339 }
340 
341 // Nonmember constructors.
342 // ------------------------------------------------------------------------
343 
344 template <class Scalar>
347 {
348  auto ter = Teuchos::rcp(new TimeEventRange<Scalar>());
349  if (pl == Teuchos::null) return ter; // Return default TimeEventRange.
350 
351  TEUCHOS_TEST_FOR_EXCEPTION(pl->get<std::string>("Type", "Range") != "Range",
352  std::logic_error,
353  "Error - Time Event Type != 'Range'. (='" +
354  pl->get<std::string>("Type") + "')\n");
355 
356  auto validPL = *ter->getValidParameters();
357  bool numEventsFound = pl->isParameter("Number of Events");
358  if (!numEventsFound) validPL.remove("Number of Events");
359 
361 
362  ter->setName(pl->get("Name", "From createTimeEventRange"));
363  ter->setTimeStart(pl->get("Start Time", ter->getTimeStart()));
364  ter->setTimeStop(pl->get("Stop Time", ter->getTimeStop()));
365  ter->setTimeStride(pl->get("Stride Time", ter->getTimeStride()));
366  if (numEventsFound)
367  ter->setNumEvents(pl->get("Number of Events", ter->getNumEvents()));
368  ter->setRelTol(pl->get("Relative Tolerance", ter->getRelTol()));
369  ter->setLandOnExactly(pl->get("Land On Exactly", ter->getLandOnExactly()));
370 
371  return ter;
372 }
373 
374 } // namespace Tempus
375 #endif // Tempus_TimeEventRange_impl_hpp
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Return a valid ParameterList with current settings.
virtual void setType(std::string s)
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
Describe member data.
virtual void setRelTol(Scalar relTol)
Set the relative tolerance.
Scalar stride_
Stride of time range.
virtual void setNumEvents(int numEvents)
Set the number of time events.
virtual void setName(std::string name)
Set the name of the TimeEvent.
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)
virtual void setLandOnExactly(bool LOE)
Set if the time event should be landed on exactly.
Scalar start_
Start of time range.
bool isParameter(const std::string &name) const
virtual bool eventInRange(Scalar time1, Scalar time2) const
Test if an event occurs within the time range.
virtual Scalar getDefaultTol() const
Return the default tolerance used by TimeEvents.
TimeEventRange specifies a start, stop and stride time.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual Scalar timeToNextEvent(Scalar time) const
How much time until the next event.
Teuchos::RCP< TimeEventRange< Scalar > > createTimeEventRange(Teuchos::RCP< Teuchos::ParameterList > pList)
Nonmember Constructor via ParameterList.
virtual void setTimeStop(Scalar stop)
Set the stop of the time range.
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
virtual void setTimeScale()
Set the time scale for the time events.
static magnitudeType magnitude(T a)
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
bool approxEqualAbsTol(Scalar a, Scalar b, Scalar absTol)
Test if values are approximately equal within the absolute tolerance.
bool approxZero(Scalar value, Scalar tol=Teuchos::ScalarTraits< Scalar >::sfmin())
Test if value is approximately zero within tolerance.
virtual void setTimeStart(Scalar start)
Set the start of the time range.
Scalar stop_
Stop of time range.
ParameterList & setName(const std::string &name)
virtual void setTimeStride(Scalar stride)
Set the stride of the time range.
virtual bool isTime(Scalar time) const
Test if time is near an event (within tolerance).
TimeEventRange()
Default constructor.
virtual Scalar timeOfNextEvent(Scalar time) const
Return the time of the next event following the input time.
virtual void setTimeRange(Scalar start, Scalar stop, Scalar stride)
Set the range of time events from start, stop and stride.