Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Tempus_SolutionHistory_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_SolutionHistory_impl_hpp
10 #define Tempus_SolutionHistory_impl_hpp
11 
12 // Teuchos
13 #include "Teuchos_StandardParameterEntryValidators.hpp"
14 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
15 #include "Teuchos_TimeMonitor.hpp"
16 
17 // Tempus
18 #include "Tempus_SolutionStateMetaData.hpp"
20 
21 //#include "Thyra_VectorStdOps.hpp"
22 
23 
24 namespace {
25 
26  static std::string Invalid_name = "Invalid";
27  static std::string KeepNewest_name = "Keep Newest";
28  static std::string Undo_name = "Undo";
29  static std::string Static_name = "Static";
30  static std::string Unlimited_name = "Unlimited";
31  static std::string Storage_name = "Storage Type";
32  static std::string Storage_default = Undo_name;
33 
34  static std::string StorageLimit_name = "Storage Limit";
35  static int StorageLimit_default = 2;
36 
37  std::vector<std::string> HistoryPolicies =
38  {Invalid_name, KeepNewest_name, Undo_name, Static_name, Unlimited_name};
39 
40  const Teuchos::RCP<Teuchos::StringToIntegralParameterEntryValidator<Tempus::StorageType> >
41  StorageTypeValidator = Teuchos::rcp(
42  new Teuchos::StringToIntegralParameterEntryValidator<Tempus::StorageType>(
43  HistoryPolicies,
44  Teuchos::tuple<Tempus::StorageType>(
50  Storage_name));
51 
52 } // namespace
53 
54 
55 namespace Tempus {
56 
57 template<class Scalar>
59  Teuchos::RCP<Teuchos::ParameterList> pList)
60 {
61  using Teuchos::RCP;
62  // Create history, an array of solution states.
63  history_ = rcp(new std::vector<RCP<SolutionState<Scalar> > >);
64 
65  this->setParameterList(pList);
66 
67  if (Teuchos::as<int>(this->getVerbLevel()) >=
68  Teuchos::as<int>(Teuchos::VERB_HIGH)) {
69  RCP<Teuchos::FancyOStream> out = this->getOStream();
70  Teuchos::OSTab ostab(out,1,"SolutionHistory::SolutionHistory");
71  *out << this->description() << std::endl;
72  }
73 }
74 
75 
76 template<class Scalar>
78  const Teuchos::RCP<SolutionState<Scalar> >& state)
79 {
80  // Check that we're not going to exceed our storage limit:
81  if (Teuchos::as<int>(history_->size()+1) > storageLimit_) {
82  switch (storageType_) {
83  case STORAGE_TYPE_INVALID: {
84  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
85  "Error - Storage type is STORAGE_TYPE_INVALID.\n");
86  break;
87  }
90  case STORAGE_TYPE_UNDO: {
91  if (state->getTime() >= history_->front()->getTime()) {
92  // Case: State is older than the youngest state in history.
93  // Remove state from the beginning of history, then add new state.
94  history_->erase(history_->begin());
95  } else {
96  // Case: State is younger than the youngest state in history.
97  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
98  Teuchos::OSTab ostab(out,1,"SolutionHistory::addState");
99  *out << "Warning, state is younger than youngest state in history. "
100  << "State not added!" << std::endl;
101  return;
102  }
103  break;
104  }
106  break;
107  default:
108  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
109  "Error - unknown storage type.\n");
110  }
111  }
112 
113  // Add new state in chronological order.
114  if (history_->size() == 0) {
115  history_->push_back(state);
116  } else {
117  typename std::vector<Teuchos::RCP<SolutionState<Scalar> > >::iterator
118  state_it = history_->begin();
119  for (; state_it < history_->end(); state_it++) {
120  if (state->getTime() < (*state_it)->getTime()) break;
121  }
122  history_->insert(state_it, state);
123  }
124 
125  TEUCHOS_TEST_FOR_EXCEPTION(getNumStates() <= 0, std::logic_error,
126  "Error - SolutionHistory::addState() Invalid history size!\n");
127 
128  return;
129 }
130 
131 template<class Scalar>
133  const Teuchos::RCP<SolutionState<Scalar> >& state, const bool updateTime)
134 {
135  using Teuchos::RCP;
136 
137  addState(state);
138  workingState_ = (*history_)[getNumStates()-1];
139  RCP<SolutionStateMetaData<Scalar> > csmd = getCurrentState()->getMetaData();
140  RCP<SolutionStateMetaData<Scalar> > wsmd = workingState_ ->getMetaData();
141  wsmd->setSolutionStatus(Status::WORKING);
142  wsmd->setIStep(csmd->getIStep()+1);
143  if (updateTime) {
144  wsmd->setTime(csmd->getTime() + csmd->getDt());
145  wsmd->setDt(csmd->getDt());
146  }
147 }
148 
149 template<class Scalar>
151  const Teuchos::RCP<SolutionState<Scalar> >& state)
152 {
153  if (history_->size() != 0) {
154  auto state_it = history_->rbegin();
155  for ( ; state_it < history_->rend(); state_it++) {
156  if (state->getTime() == (*state_it)->getTime()) break;
157  }
158 
159  TEUCHOS_TEST_FOR_EXCEPTION(state_it == history_->rend(), std::logic_error,
160  "Error - removeState() Could not remove state = "
161  // << state_it->describe()
162  );
163 
164  // Need to be careful when erasing a reverse iterator.
165  history_->erase(std::next(state_it).base());
166  }
167  return;
168 }
169 
170 
171 template<class Scalar>
173 {
174  Teuchos::RCP<SolutionState<Scalar> > tmpState = findState(time);
175  removeState(tmpState);
176 }
177 
178 
179 template<class Scalar>
180 Teuchos::RCP<SolutionState<Scalar> >
181 SolutionHistory<Scalar>::findState(const Scalar time) const
182 {
183  TEUCHOS_TEST_FOR_EXCEPTION(
184  !(minTime() <= time and time <= maxTime()), std::logic_error,
185  "Error - SolutionHistory::findState() Requested time out of range!\n"
186  " [Min, Max] = [" << minTime() << ", " << maxTime() << "]\n"
187  " time = "<< time <<"\n");
188 
189  // Use last step in solution history as the scale for comparing times
190  const Scalar scale =
191  history_->size() > 0 ? (*history_)[history_->size()-1]->getTime() : Scalar(1.0);
192  // Linear search
193  auto state_it = history_->begin();
194  for ( ; state_it < history_->end(); ++state_it) {
195  if (floating_compare_equals((*state_it)->getTime(),time,scale))
196  break;
197  }
198 
199  TEUCHOS_TEST_FOR_EXCEPTION(state_it == history_->end(), std::logic_error,
200  "Error - SolutionHistory::findState()!\n"
201  " Did not find a SolutionState with time = " <<time<< std::endl);
202 
203  return *state_it;
204 }
205 
206 
207 template<class Scalar>
208 Teuchos::RCP<SolutionState<Scalar> >
210 {
211  Teuchos::RCP<SolutionState<Scalar> > state_out = getCurrentState()->clone();
212  interpolate<Scalar>(*interpolator_, history_, time, state_out.get());
213  return state_out;
214 }
215 
216 
217 template<class Scalar>
218 void
220  const Scalar time, SolutionState<Scalar>* state_out) const
221 {
222  interpolate<Scalar>(*interpolator_, history_, time, state_out);
223 }
224 
225 
226 /** Initialize the working state */
227 template<class Scalar>
229 {
230  TEMPUS_FUNC_TIME_MONITOR("Tempus::SolutionHistory::initWorkingState()");
231  {
232  TEUCHOS_TEST_FOR_EXCEPTION(getCurrentState() == Teuchos::null,
233  std::logic_error,
234  "Error - SolutionHistory::initWorkingState()\n"
235  "Can not initialize working state without a current state!\n");
236 
237  // If workingState_ has a valid pointer, we are still working on it,
238  // i.e., step failed and trying again, so do not initialize it.
239  if (getWorkingState(false) != Teuchos::null) return;
240 
241  Teuchos::RCP<SolutionState<Scalar> > newState;
242  if (getNumStates() < storageLimit_) {
243  // Create newState which is duplicate of currentState
244  newState = getCurrentState()->clone();
245  } else {
246  // Recycle old state and copy currentState
247  newState = (*history_)[0];
248  history_->erase(history_->begin());
249  if (getNumStates() > 0) newState->copy(getCurrentState());
250  // When using the Griewank algorithm, we will want to select which
251  // older state to recycle.
252  }
253 
254  addWorkingState(newState);
255 
256  }
257  return;
258 }
259 
260 
261 template<class Scalar>
263 {
264  Teuchos::RCP<SolutionStateMetaData<Scalar> > md =
265  getWorkingState()->getMetaData();
266 
267  if ( md->getSolutionStatus() == Status::PASSED ) {
268  md->setNFailures(std::max(0,md->getNFailures()-1));
269  md->setNConsecutiveFailures(0);
270  md->setSolutionStatus(Status::PASSED);
271  //md->setIsSynced(true);
272  md->setIsInterpolated(false);
273  workingState_ = Teuchos::null;
274  } else {
275  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
276  Teuchos::OSTab ostab(out,1,"SolutionHistory::promoteWorkingState()");
277  *out << "Warning - WorkingState is not passing, so not promoted!\n"
278  << std::endl;
279  }
280 }
281 
282 
283 template<class Scalar>
285 {
286  storageLimit_ = std::max(1,storage_limit);
287 
288  TEUCHOS_TEST_FOR_EXCEPTION(
289  (Teuchos::as<int>(history_->size()) > storageLimit_), std::logic_error,
290  "Error - requested storage limit = " << storageLimit_
291  << " is smaller than the current number of states stored = "
292  << history_->size() << "!\n");
293 }
294 
295 
296 template<class Scalar>
297 Teuchos::RCP<SolutionState<Scalar> >
299 {
300  const int m = history_->size();
301  TEUCHOS_TEST_FOR_EXCEPTION( (m < 1), std::out_of_range,
302  "Error - getStateTimeIndexN() No states in SolutionHistory!\n");
303  return (*history_)[m-1];
304 }
305 
306 
307 template<class Scalar>
308 Teuchos::RCP<SolutionState<Scalar> >
310 {
311  const int m = history_->size();
312  TEUCHOS_TEST_FOR_EXCEPTION( (m < 2), std::out_of_range,
313  "Error - getStateTimeIndexNM1() Not enough states in "
314  << "SolutionHistory!\n");
315  const int n = (*history_)[m-1]->getIndex();
316  const int nm1 = (*history_)[m-2]->getIndex();
317 
318  // No need to search SolutionHistory as states n and nm1 should be
319  // next to each other.
320  TEUCHOS_TEST_FOR_EXCEPTION( (nm1 != n-1), std::out_of_range,
321  "Error - getStateTimeIndexNM1() Timestep index n-1 is not in "
322  << "SolutionHistory!\n"
323  << " (n)th index = " << n << "\n"
324  << " (n-1)th index = " << nm1 << "\n");
325 
326  return (*history_)[m-2];
327 }
328 
329 
330 template<class Scalar>
331 Teuchos::RCP<SolutionState<Scalar> >
333 {
334  const int m = history_->size();
335  TEUCHOS_TEST_FOR_EXCEPTION( (m < 3), std::out_of_range,
336  "Error - getStateTimeIndexNM1() Not enough states in "
337  << "SolutionHistory!\n");
338  const int n = (*history_)[m-1]->getIndex();
339  const int nm2 = (*history_)[m-3]->getIndex();
340 
341  // Assume states n and nm2 are one away from each other.
342  // May need to do a search otherwise.
343  TEUCHOS_TEST_FOR_EXCEPTION( (nm2 != n-2), std::out_of_range,
344  "Error - getStateTimeIndexNM1() Timestep index n-2 is not in "
345  << "SolutionHistory!\n"
346  << " (n)th index = " << n << "\n"
347  << " (n-2)th index = " << nm2 << "\n");
348 
349  return (*history_)[m-3];
350 }
351 
352 
353 template<class Scalar>
354 Teuchos::RCP<SolutionState<Scalar> >
356 {
357  typename std::vector<Teuchos::RCP<SolutionState<Scalar> > >::iterator
358  state_it = history_->begin();
359  for (; state_it < history_->end(); state_it++) {
360  if ((*state_it)->getIndex() == index) break;
361  }
362  TEUCHOS_TEST_FOR_EXCEPTION( state_it==history_->end(), std::out_of_range,
363  "Error - getStateTimeIndex() Timestep index is not in "
364  << "SolutionHistory!\n"
365  << " index = " << index << "\n");
366  return (*state_it);
367 }
368 
369 
370 template<class Scalar>
372 {
373  return ("Tempus::SolutionHistory - name = '" + name_ + "'");
374 }
375 
376 
377 template<class Scalar>
379  Teuchos::FancyOStream &out,
380  const Teuchos::EVerbosityLevel verbLevel) const
381 {
382  if ((Teuchos::as<int>(verbLevel)==Teuchos::as<int>(Teuchos::VERB_DEFAULT)) ||
383  (Teuchos::as<int>(verbLevel)>=Teuchos::as<int>(Teuchos::VERB_LOW) ) ){
384  out << description() << "::describe" << std::endl;
385  //out << "interpolator = " << interpolator->description() << std::endl;
386  out << "storageLimit = " << storageLimit_ << std::endl;
387  out << "storageType = " << storageType_ << std::endl;
388  out << "number of states = " << history_->size() << std::endl;
389  out << "time range = (" << history_->front()->getTime() << ", "
390  << history_->back()->getTime() << ")"
391  << std::endl;
392  } else if (Teuchos::as<int>(verbLevel) >=
393  Teuchos::as<int>(Teuchos::VERB_HIGH)) {
394  out << "SolutionStates: " << std::endl;
395  for (int i=0; i<(int)history_->size() ; ++i) {
396  out << "SolutionState[" << i << "] = " << std::endl;
397  (*history_)[i]->describe(out,this->getVerbLevel());
398  }
399  }
400 }
401 
402 
403 template <class Scalar>
405  Teuchos::RCP<Teuchos::ParameterList> const& pList)
406 {
407  if (pList == Teuchos::null) {
408  // Create default parameters if null, otherwise keep current parameters.
409  if (shPL_ == Teuchos::null) {
410  shPL_ = Teuchos::parameterList("Solution History");
411  *shPL_ = *(this->getValidParameters());
412  }
413  } else {
414  shPL_ = pList;
415  }
416  shPL_->validateParametersAndSetDefaults(*this->getValidParameters());
417 
418  name_ = shPL_->name();
419 
420  storageType_ = StorageTypeValidator->getIntegralValue(
421  *shPL_, Storage_name, Storage_default);
422 
423  int storage_limit = shPL_->get(StorageLimit_name, StorageLimit_default);
424 
425  switch (storageType_) {
428  storageType_ = STORAGE_TYPE_KEEP_NEWEST;
429  if (storage_limit != 1) {
430  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
431  Teuchos::OSTab ostab(out,1,"SolutionHistory::setParameterList");
432  *out << "Warning - 'Storage Limit' for 'Keep Newest' is 1.\n"
433  << " (Storage Limit = "<<storage_limit<<"). Resetting to 1."
434  << std::endl;
435  storage_limit = 1;
436  }
437  setStorageLimit(storage_limit);
438  break;
439  }
440  case STORAGE_TYPE_UNDO: {
441  if (storage_limit != 2) {
442  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
443  Teuchos::OSTab ostab(out,1,"SolutionHistory::setParameterList");
444  *out << "Warning - 'Storage Limit' for 'Undo' is 2.\n"
445  << " (Storage Limit = "<<storage_limit<<"). Resetting to 2."
446  << std::endl;
447  storage_limit = 2;
448  }
449  setStorageLimit(storage_limit);
450  break;
451  }
452  case STORAGE_TYPE_STATIC: {
453  break;
454  }
455  case STORAGE_TYPE_UNLIMITED: {
456  storage_limit = 1000000000;
457  break;
458  }
459  }
460  setStorageLimit(storage_limit);
461 
463  Teuchos::sublist(shPL_, "Interpolator"));
464 }
465 
466 
467 template<class Scalar>
468 Teuchos::RCP<const Teuchos::ParameterList>
470 {
471  Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
472 
473  pl->setName("Valid ParameterList");
474 
475  pl->set(Storage_name, Storage_default,
476  "'Storage Type' sets the memory storage. "
477  "'Keep Newest' - will retain the single newest solution state. "
478  "'Undo' - will retain two solution states in order to do a single undo. "
479  "'Static' - will retain 'Storage Limit' number of solution states. "
480  "'Unlimited' - will not remove any solution states!",
481  StorageTypeValidator);
482 
483  pl->set(StorageLimit_name, StorageLimit_default,
484  "Storage limit for the solution history.");
485 
486  // Interpolator
487  pl->sublist("Interpolator",false,"").disableRecursiveValidation();
488 
489  return pl;
490 }
491 
492 
493 template <class Scalar>
494 Teuchos::RCP<Teuchos::ParameterList>
496 {
497  return(shPL_);
498 }
499 
500 
501 template <class Scalar>
502 Teuchos::RCP<Teuchos::ParameterList>
504 {
505  Teuchos::RCP<Teuchos::ParameterList> temp_plist = shPL_;
506  shPL_ = Teuchos::null;
507  return(temp_plist);
508 }
509 
510 // Nonmember constructor.
511 template<class Scalar>
512 Teuchos::RCP<SolutionHistory<Scalar> > solutionHistory(
513  Teuchos::RCP<Teuchos::ParameterList> pList)
514 {
515  Teuchos::RCP<SolutionHistory<Scalar> > sh=rcp(new SolutionHistory<Scalar>(pList));
516  return sh;
517 }
518 
519 template<class Scalar>
521  const Teuchos::RCP<Interpolator<Scalar> >& interpolator)
522 {
523  if (interpolator == Teuchos::null) {
525  } else {
526  interpolator_ = interpolator;
527  }
528  if (Teuchos::as<int>(this->getVerbLevel()) >=
529  Teuchos::as<int>(Teuchos::VERB_HIGH)) {
530  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
531  Teuchos::OSTab ostab(out,1,"SolutionHistory::setInterpolator");
532  *out << "interpolator = " << interpolator_->description() << std::endl;
533  }
534 }
535 
536 template<class Scalar>
537 Teuchos::RCP<Interpolator<Scalar> >
539 {
540  return interpolator_;
541 }
542 
543 template<class Scalar>
544 Teuchos::RCP<const Interpolator<Scalar> >
546 {
547  return interpolator_;
548 }
549 
550 template<class Scalar>
551 Teuchos::RCP<Interpolator<Scalar> >
553 {
554  Teuchos::RCP<Interpolator<Scalar> > old_interpolator = interpolator_;
555  interpolator_ = lagrangeInterpolator<Scalar>();
556  return old_interpolator;
557 }
558 
559 
560 } // namespace Tempus
561 #endif // Tempus_SolutionHistory_impl_hpp
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndexNM2() const
Get the state with timestep index equal to n-2.
Teuchos::RCP< Interpolator< Scalar > > unSetInterpolator()
Unset the interpolator for this history.
Keep the 2 newest states for undo.
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void initWorkingState()
Initialize the working state.
Teuchos::RCP< SolutionState< Scalar > > interpolateState(const Scalar time) const
Generate and interpolate a new solution state at requested time.
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
Teuchos::RCP< Interpolator< Scalar > > getNonconstInterpolator()
void setStorageLimit(int storage_limit)
Set the maximum storage of this history.
static Teuchos::RCP< Interpolator< Scalar > > createInterpolator(std::string interpolatorType="")
Create default interpolator from interpolator type (e.g., &quot;Linear&quot;).
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndexNM1() const
Get the state with timestep index equal to n-1.
void promoteWorkingState()
Promote the working state to current state.
void removeState(const Teuchos::RCP< SolutionState< Scalar > > &state)
Remove solution state.
Base strategy class for interpolation functionality.
virtual std::string description() const
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &pl)
Teuchos::RCP< SolutionHistory< Scalar > > solutionHistory(Teuchos::RCP< Teuchos::ParameterList > pList=Teuchos::null)
Nonmember constructor.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
SolutionHistory is basically a container of SolutionStates. SolutionHistory maintains a collection of...
void addState(const Teuchos::RCP< SolutionState< Scalar > > &state)
Add solution state to history.
void addWorkingState(const Teuchos::RCP< SolutionState< Scalar > > &state, const bool updateTime=true)
Add a working solution state to history.
Keep a fix number of states.
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndex(int index) const
Get the state with timestep index equal to &quot;index&quot;.
SolutionHistory(Teuchos::RCP< Teuchos::ParameterList > shPL=Teuchos::null)
Contructor.
Teuchos::RCP< SolutionState< Scalar > > getStateTimeIndexN() const
Get the state with timestep index equal to n.
Teuchos::RCP< const Interpolator< Scalar > > getInterpolator() const
bool floating_compare_equals(const Scalar &a, const Scalar &b, const Scalar &scale)
Helper function for comparing times.
void setInterpolator(const Teuchos::RCP< Interpolator< Scalar > > &interpolator)
Set the interpolator for this history.
Teuchos::RCP< SolutionState< Scalar > > findState(const Scalar time) const
Find solution state at requested time (no interpolation)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...