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()),
29  setTimeRange(0.0, 0.0, 0.0);
30  setLandOnExactly(true);
31  std::ostringstream oss;
32  oss << "TimeEventRange (" << start_<< "; " << stop_<< "; " << stride_<< ")";
33  this->setName(oss.str());
34 }
35 
36 
37 template<class Scalar>
39  Scalar start, Scalar stop, Scalar stride,
40  std::string name, bool landOnExactly, Scalar relTol)
41  : start_ (start),
42  stop_ (stop),
43  stride_(stride),
44  numEvents_(0),
45  timeScale_(std::max(std::abs(start_), std::abs(stop_))),
46  relTol_(relTol),
47  absTol_(relTol_*timeScale_),
48  landOnExactly_(landOnExactly)
49 {
50  this->setType("Range");
51  if (name == "") {
52  std::ostringstream oss;
53  oss << "TimeEventRange (" << start << "; " << stop << "; " << stride << ")";
54  this->setName(oss.str());
55  } else {
56  this->setName(name);
57  }
58  setRelTol(relTol);
59  setTimeRange(start, stop, stride);
60  setLandOnExactly(landOnExactly);
61 }
62 
63 
64 template<class Scalar>
66  Scalar start, Scalar stop, int numEvents,
67  std::string name, bool landOnExactly, Scalar relTol)
68  : start_ (start),
69  stop_ (stop),
70  stride_(0.0),
71  numEvents_(numEvents),
72  timeScale_(std::max(std::abs(start_), std::abs(stop_))),
73  relTol_(relTol),
74  absTol_(relTol_*timeScale_),
75  landOnExactly_(landOnExactly)
76 {
77  if (name == "") {
78  std::ostringstream oss;
79  oss << "TimeEventRange (" << start << "; " << stop << "; " << numEvents << ")";
80  this->setName(oss.str());
81  } else {
82  this->setName(name);
83  }
84  setRelTol(relTol);
85  setTimeRange(start, stop, numEvents);
86  setLandOnExactly(landOnExactly);
87 }
88 
89 
90 template<class Scalar>
92  Scalar start, Scalar stop, 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 
106 template<class Scalar>
108  Scalar start, Scalar stop, int numEvents)
109 {
110  start_ = start;
111  stop_ = stop;
112  if (stop_ < start_) {
113  Scalar tmp = start_;
114  start_ = stop_;
115  stop_ = tmp;
116  }
117  setTimeScale();
118  setNumEvents(numEvents);
119 }
120 
121 
122 template<class Scalar>
124 {
125  start_ = start;
126  if (stop_ < start_) stop_ = start_;
127  setTimeScale();
128  setTimeStride(stride_); // Reset numEvents with the current stride.
129 }
130 
131 
132 template<class Scalar>
134 {
135  stop_ = stop;
136  if (start_ > stop_) start_ = stop_;
137  setTimeScale();
138  setTimeStride(stride_); // Reset numEvents with the current stride.
139 }
140 
141 
142 template<class Scalar>
144 {
145  timeScale_ = std::max(std::abs(start_), std::abs(stop_));
146  absTol_ = relTol_*timeScale_;
147 
148  // Check if timeScale is near zero.
149  if ( approxZero(timeScale_, absTol_) ) {
150  timeScale_ = 1.0;
151  absTol_ = relTol_*timeScale_;
152  }
153 }
154 
155 
156 template<class Scalar>
158 {
160  if ( approxEqualAbsTol(start_, stop_, absTol_) ) {
161  stride_ = 0.0;
162  numEvents_ = 1;
163  return;
164  }
165 
166  if ((stride_ > stop_ - start_) || (stride_ < 2*absTol_)) {
167  stride_ = stop_ - start_;
168  }
169 
170  numEvents_ = int((stop_+absTol_ - start_) / stride_) + 1;
171 }
172 
173 
174 template<class Scalar>
176 {
177  numEvents_ = numEvents;
178  if (numEvents_ < 0) { // Do not use numEvents_ to set stride! Reset numEvents_ from stride_.
179  if (stride_ < 2 * absTol_) stride_ = 2*absTol_;
180  numEvents_ = int((stop_+absTol_ - start_) / stride_) + 1;
181  return;
182  } else if ( approxEqualAbsTol(start_, stop_, absTol_) ) {
183  numEvents_ = 1;
184  stride_ = 0.0;
185  stride_ = stop_ - start_;
186  } else {
187  if (numEvents_ < 2) numEvents_ = 2;
188  stride_ = (stop_ - start_)/Scalar(numEvents_-1);
189  }
190 
191  // If stride_ is smaller than twice the absolute tolerance,
192  // the time steps cannot be differentiated!
193  if (stride_ <= 2 * absTol_) setTimeStride(2*absTol_);
194 }
195 
196 
197 template<class Scalar>
199 {
200  relTol_ = std::abs(relTol);
201  setTimeScale();
202 }
203 
204 
205 template<class Scalar>
206 bool TimeEventRange<Scalar>::isTime(Scalar time) const
207 {
208  // Check if before first event.
209  if (time < start_-absTol_) return false;
210 
211  // Check if after last event.
212  const Scalar timeOfLast = start_ + (numEvents_-1) * stride_;
213  if (time > timeOfLast+absTol_) return false;
214 
215  int numStrides = 0;
216  if ( !approxZero(stride_, 2*absTol_) )
217  numStrides = (time - start_) / stride_;
218 
219  numStrides = std::min(std::max(0, numStrides), int(numEvents_-1));
220  const Scalar leftBracket = start_ + numStrides * stride_;
221 
222  // Check if close to left bracket.
223  if ( approxEqualAbsTol(time, leftBracket, absTol_) )
224  return true;
225 
226  // Check if close to right bracket.
227  const Scalar rightBracket = leftBracket + stride_;
228  if ( approxEqualAbsTol(time, rightBracket, absTol_) )
229  return true;
230 
231  return false;
232 }
233 
234 
235 template<class Scalar>
237 {
238  return timeOfNextEvent(time) - time; // Neg. time indicates in the past.
239 }
240 
241 
242 template<class Scalar>
244 {
245  // Check if before first event.
246  if (time < start_-absTol_) return start_;
247 
248  const Scalar timeOfLast = start_ + (numEvents_-1) * stride_;
249  // Check if after or close to last event.
250  if (time > timeOfLast-absTol_) return std::numeric_limits<Scalar>::max();
251 
252  int numStrides = 0;
253  if ( !approxZero(stride_, 2*absTol_) )
254  numStrides = int((time - start_) / stride_) + 1;
255  const Scalar timeEvent = start_ + numStrides*stride_;
256 
257  // Check timeEvent is near time. If so, return the next event.
258  if ( approxEqualAbsTol(time, timeEvent, absTol_) )
259  return timeEvent + stride_;
260 
261  return timeEvent;
262 }
263 
264 
265 template<class Scalar>
266 bool TimeEventRange<Scalar>::eventInRange(Scalar time1, Scalar time2) const
267 {
268  if (time1 > time2) {
269  Scalar tmp = time1;
270  time1 = time2;
271  time2 = tmp;
272  }
273 
274  // Check if range is completely outside time events.
275  const Scalar timeOfLast = start_ + (numEvents_-1) * stride_;
276  if (time2 < start_-absTol_ || timeOfLast+absTol_ < time1) return false;
277 
278  if ( approxZero(stride_) )
279  return (time1 < start_-absTol_ && start_-absTol_ <= time2);
280 
281  const int strideJustBeforeTime1 = std::min(int(numEvents_-1),
282  std::max(int(0), int((time1 - start_ + absTol_) / stride_ - 0.5)));
283 
284  const int strideJustAfterTime2 = std::max(int(0), std::min(int(numEvents_-1),
285  int((time2 - start_ + absTol_) / stride_ + 0.5)));
286 
287  for ( int i = strideJustBeforeTime1; i <= strideJustAfterTime2; i++ ) {
288  const Scalar timeEvent = start_ + i * stride_;
289  if (time1 < timeEvent-absTol_ && timeEvent-absTol_ <= time2) return true;
290  }
291 
292  return false;
293 }
294 
295 
296 template<class Scalar>
298  const Teuchos::EVerbosityLevel verbLevel) const
299 {
300  auto l_out = Teuchos::fancyOStream( out.getOStream() );
301  Teuchos::OSTab ostab(*l_out, 2, "TimeEventRange");
302  l_out->setOutputToRootOnly(0);
303 
304  *l_out << "TimeEventRange:" << "\n"
305  << " name = " << this->getName() << "\n"
306  << " Type = " << this->getType() << "\n"
307  << " start_ = " << start_ << "\n"
308  << " stop_ = " << stop_ << "\n"
309  << " stride_ = " << stride_ << "\n"
310  << " numEvents_ = " << numEvents_ << "\n"
311  << " timeScale_ = " << timeScale_ << "\n"
312  << " relTol_ = " << relTol_ << "\n"
313  << " absTol_ = " << absTol_ << "\n"
314  << " landOnExactly_ = " << landOnExactly_ << std::endl;
315 }
316 
317 
318 template<class Scalar>
321 {
323  Teuchos::parameterList("Time Event Range");
324 
325  pl->setName(this->getName());
326  pl->set("Name", this->getName());
327  pl->set("Type", this->getType());
328 
329  pl->set("Start Time", getTimeStart(), "Start of time range");
330  pl->set("Stop Time", getTimeStop(), "Stop of time range");
331  pl->set("Stride Time", getTimeStride(), "Stride of time range");
332 
333  if ( getTimeStride()*Scalar(getNumEvents()-1) - (getTimeStop()-getTimeStart()) < getAbsTol() )
334  pl->set("Number of Events", getNumEvents(),
335  "Number of events in time range. If specified, 'Stride Time' is reset.");
336 
337  pl->set("Relative Tolerance", getRelTol(),
338  "Relative time tolerance for matching time events.");
339 
340  pl->set("Land On Exactly", getLandOnExactly(),
341  "Should these time events be landed on exactly, i.e, adjust the timestep to hit time event, versus stepping over and keeping the time step unchanged.");
342 
343  return pl;
344 }
345 
346 
347 // Nonmember constructors.
348 // ------------------------------------------------------------------------
349 
350 template<class Scalar>
353 {
354  auto ter = Teuchos::rcp(new TimeEventRange<Scalar>());
355  if (pl == Teuchos::null) return ter; // Return default TimeEventRange.
356 
358  pl->get<std::string>("Type", "Range") != "Range",
359  std::logic_error,
360  "Error - Time Event Type != 'Range'. (='"
361  + pl->get<std::string>("Type")+"')\n");
362 
363  auto validPL = *ter->getValidParameters();
364  bool numEventsFound = pl->isParameter("Number of Events");
365  if (!numEventsFound) validPL.remove("Number of Events");
366 
368 
369  ter->setName (pl->get("Name", "From createTimeEventRange"));
370  ter->setTimeStart (pl->get("Start Time", ter->getTimeStart()));
371  ter->setTimeStop (pl->get("Stop Time", ter->getTimeStop()));
372  ter->setTimeStride (pl->get("Stride Time", ter->getTimeStride()));
373  if (numEventsFound)
374  ter->setNumEvents (pl->get("Number of Events", ter->getNumEvents()));
375  ter->setRelTol (pl->get("Relative Tolerance", ter->getRelTol()));
376  ter->setLandOnExactly (pl->get("Land On Exactly", ter->getLandOnExactly()));
377 
378  return ter;
379 }
380 
381 
382 } // namespace Tempus
383 #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.