Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_SolutionState_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_SolutionState_impl_hpp
10 #define Tempus_SolutionState_impl_hpp
11 
12 #include "Thyra_VectorStdOps.hpp"
13 
14 namespace Tempus {
15 
16 
17 template<class Scalar>
19  : x_ (Teuchos::null),
20  x_nc_ (Teuchos::null),
21  xdot_ (Teuchos::null),
22  xdot_nc_ (Teuchos::null),
23  xdotdot_ (Teuchos::null),
24  xdotdot_nc_ (Teuchos::null),
25  stepperState_ (Teuchos::null),
26  stepperState_nc_(Teuchos::null),
27  physicsState_ (Teuchos::null),
28  physicsState_nc_(Teuchos::null)
29 {
36 }
37 
38 
39 #ifndef TEMPUS_HIDE_DEPRECATED_CODE
40 template<class Scalar>
44  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot,
45  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
46  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
47  : x_ (x),
48  x_nc_ (x),
49  xdot_ (xdot),
50  xdot_nc_ (xdot),
51  xdotdot_ (xdotdot),
52  xdotdot_nc_ (xdotdot),
53  stepperState_ (stepperState),
54  stepperState_nc_(stepperState),
55  physicsState_ (physicsState),
56  physicsState_nc_(physicsState)
57 {
60  if (stepperState_nc_ == Teuchos::null) {
63  }
64  if (physicsState_nc_ == Teuchos::null) {
67  }
68 }
69 
70 template<class Scalar>
72  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
73  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
74  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot,
75  const Teuchos::RCP<const StepperState<Scalar> >& stepperState,
76  const Teuchos::RCP<const PhysicsState<Scalar> >& physicsState)
77  : x_ (x),
78  x_nc_ (Teuchos::null),
79  xdot_ (xdot),
80  xdot_nc_ (Teuchos::null),
81  xdotdot_ (xdotdot),
82  xdotdot_nc_ (Teuchos::null),
83  stepperState_ (stepperState),
84  stepperState_nc_(Teuchos::null),
85  physicsState_ (physicsState),
86  physicsState_nc_(Teuchos::null)
87 {
90 
91  using Teuchos::rcp_const_cast;
92  if (stepperState_ == Teuchos::null) {
94  }
95  if (physicsState_ == Teuchos::null) {
97  }
98 }
99 #endif
100 
101 
102 template<class Scalar>
107  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot,
108  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
109  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
110  : metaData_ (metaData),
111  metaData_nc_ (metaData),
112  x_ (x),
113  x_nc_ (x),
114  xdot_ (xdot),
115  xdot_nc_ (xdot),
116  xdotdot_ (xdotdot),
117  xdotdot_nc_ (xdotdot),
118  stepperState_ (stepperState),
119  stepperState_nc_(stepperState),
120  physicsState_ (physicsState),
121  physicsState_nc_(physicsState)
122 {
123  if (stepperState_nc_ == Teuchos::null) {
126  }
127  if (physicsState_nc_ == Teuchos::null) {
130  }
131 }
132 
133 template<class Scalar>
135  const Teuchos::RCP<const SolutionStateMetaData<Scalar> > metaData,
136  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
137  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
138  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot,
139  const Teuchos::RCP<const StepperState<Scalar> >& stepperState,
140  const Teuchos::RCP<const PhysicsState<Scalar> >& physicsState)
141  : metaData_ (metaData),
142  metaData_nc_ (Teuchos::null),
143  x_ (x),
144  x_nc_ (Teuchos::null),
145  xdot_ (xdot),
146  xdot_nc_ (Teuchos::null),
147  xdotdot_ (xdotdot),
148  xdotdot_nc_ (Teuchos::null),
149  stepperState_ (stepperState),
150  stepperState_nc_(Teuchos::null),
151  physicsState_ (physicsState),
152  physicsState_nc_(Teuchos::null)
153 {
154  if (stepperState_ == Teuchos::null) {
156  }
157  if (physicsState_ == Teuchos::null) {
159  }
160 }
161 
162 
163 #ifndef TEMPUS_HIDE_DEPRECATED_CODE
164 template<class Scalar>
166  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
167  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
168  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
169 {
170  typedef Thyra::ModelEvaluatorBase MEB;
171  using Teuchos::rcp_const_cast;
172 
173  metaData_nc_ = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
174  metaData_nc_->setSolutionStatus(Status::PASSED);
175  metaData_ = metaData_nc_;
176 
177  MEB::InArgs<Scalar> inArgs = model->getNominalValues();
178 
179  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.supports(MEB::IN_ARG_x) == false,
180  std::logic_error,
181  model->description() << "does not support an x solution vector!");
182 
183  // The solution vector, x, is required (usually).
184  x_nc_ = rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x());
185  x_ = x_nc_;
186 
187  // The solution derivative, xdot, can be optional provided, based on
188  // application needs. Here we will base it on "supports" IN_ARG_x_dot.
189  // Depending on the stepper used, a temporary xdot vector may be created
190  // within the Stepper, but not moved to the SolutionState.
191  if (inArgs.supports(MEB::IN_ARG_x_dot)) {
192  xdot_nc_ = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
193  xdot_ = xdot_nc_;
194  } else {
195  xdot_nc_ = Teuchos::null;
196  xdot_ = xdot_nc_;
197  }
198 
199  // Similar as xdot.
200  if (inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
201  xdotdot_nc_ =
202  rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x_dot_dot());
203  xdotdot_ = xdotdot_nc_;
204  } else {
205  xdotdot_nc_ = Teuchos::null;
206  xdotdot_ = xdotdot_nc_;
207  }
208 
209  if (stepperState == Teuchos::null) {
210  stepperState_nc_ = Teuchos::rcp(new StepperState<Scalar> ()); // Use default
211  stepperState_ = stepperState_nc_;
212  } else {
213  stepperState_nc_ = stepperState;
214  stepperState_ = stepperState;
215  }
216 
217  if (physicsState == Teuchos::null) {
218  physicsState_nc_ = Teuchos::rcp(new PhysicsState<Scalar> ()); // Use default
219  physicsState_ = physicsState_nc_;
220  } else {
221  physicsState_nc_ = physicsState;
222  physicsState_ = physicsState;
223  }
224 }
225 #endif
226 
227 template<class Scalar>
229  :metaData_ (ss_.metaData_),
230  metaData_nc_ (ss_.metaData_nc_),
231  x_ (ss_.x_),
232  x_nc_ (ss_.x_nc_),
233  xdot_ (ss_.xdot_),
234  xdot_nc_ (ss_.xdot_nc_),
235  xdotdot_ (ss_.xdotdot_),
236  xdotdot_nc_ (ss_.xdotdot_nc_),
237  stepperState_ (ss_.stepperState_),
238  stepperState_nc_(ss_.stepperState_nc_),
239  physicsState_ (ss_.physicsState_),
240  physicsState_nc_(ss_.physicsState_nc_)
241 {}
242 
243 
244 template<class Scalar>
246 {
247  using Teuchos::RCP;
248 
249  RCP<SolutionStateMetaData<Scalar> > metaData_out;
250  if (!Teuchos::is_null(metaData_)) metaData_out = metaData_->clone();
251 
252  RCP<Thyra::VectorBase<Scalar> > x_out;
253  if (!Teuchos::is_null(x_)) x_out = x_->clone_v();
254 
255  RCP<Thyra::VectorBase<Scalar> > xdot_out;
256  if (!Teuchos::is_null(xdot_)) xdot_out = xdot_->clone_v();
257 
258  RCP<Thyra::VectorBase<Scalar> > xdotdot_out;
259  if (!Teuchos::is_null(xdotdot_)) xdotdot_out = xdotdot_->clone_v();
260 
261  RCP<StepperState<Scalar> > sS_out;
262  if (!Teuchos::is_null(stepperState_)) sS_out=stepperState_->clone();
263 
264  RCP<PhysicsState<Scalar> > pS_out;
265  if (!Teuchos::is_null(physicsState_)) pS_out=physicsState_->clone();
266 
267  RCP<SolutionState<Scalar> > ss_out = Teuchos::rcp(new SolutionState<Scalar> (
268  metaData_out, x_out, xdot_out, xdotdot_out, sS_out, pS_out));
269 
270  return ss_out;
271 }
272 
273 
274 template<class Scalar>
277 {
278  metaData_nc_->copy(ss->metaData_);
279  this->copySolutionData(ss);
280 }
281 
282 
283 template<class Scalar>
286 {
287  if (ss->x_ == Teuchos::null)
288  x_nc_ = Teuchos::null;
289  else {
290  if (x_nc_ == Teuchos::null) {
291  x_nc_ = ss->x_->clone_v();
292  }
293  else
294  Thyra::V_V(x_nc_.ptr(), *(ss->x_));
295  }
296  x_ = x_nc_;
297 
298  if (ss->xdot_ == Teuchos::null)
299  xdot_nc_ = Teuchos::null;
300  else {
301  if (xdot_nc_ == Teuchos::null)
302  xdot_nc_ = ss->xdot_->clone_v();
303  else
304  Thyra::V_V(xdot_nc_.ptr(), *(ss->xdot_));
305  }
306  xdot_ = xdot_nc_;
307 
308  if (ss->xdotdot_ == Teuchos::null)
309  xdotdot_nc_ = Teuchos::null;
310  else {
311  if (xdotdot_nc_ == Teuchos::null)
312  xdotdot_nc_ = ss->xdotdot_->clone_v();
313  else
314  Thyra::V_V(xdotdot_nc_.ptr(), *(ss->xdotdot_));
315  }
316  xdotdot_ = xdotdot_nc_;
317 
318  if (ss->stepperState_ == Teuchos::null)
319  stepperState_nc_ = Teuchos::null;
320  else {
321  if (stepperState_nc_ == Teuchos::null)
322  stepperState_nc_ = ss->stepperState_->clone();
323  else
324  stepperState_nc_->copy(ss->stepperState_);
325  }
326  stepperState_ = stepperState_nc_;
327 
328  if (ss->physicsState_ == Teuchos::null)
329  physicsState_nc_ = Teuchos::null;
330  else {
331  if (physicsState_nc_ == Teuchos::null)
332  physicsState_nc_ = ss->physicsState_->clone();
333  else
334  physicsState_nc_->copy(ss->physicsState_);
335  }
336  physicsState_ = physicsState_nc_;
337 }
338 
339 template<class Scalar>
341 {
342  return (this->metaData_->getTime() < ss.metaData_->getTime());
343 }
344 
345 template<class Scalar>
347 {
348  return (this->metaData_->getTime() <= ss.metaData_->getTime());
349 }
350 
351 template<class Scalar>
352 bool SolutionState<Scalar>::operator< (const Scalar& t) const
353 {
354  return (this->metaData_->getTime() < t);
355 }
356 
357 template<class Scalar>
358 bool SolutionState<Scalar>::operator<= (const Scalar& t) const
359 {
360  return (this->metaData_->getTime() <= t);
361 }
362 
363 template<class Scalar>
365 {
366  return (this->metaData_->getTime() > ss.metaData_->getTime());
367 }
368 
369 template<class Scalar>
371 {
372  return (this->metaData_->getTime() >= ss.metaData_->getTime());
373 }
374 
375 template<class Scalar>
376 bool SolutionState<Scalar>::operator> (const Scalar& t) const
377 {
378  return (this->metaData_->getTime() > t);
379 }
380 
381 template<class Scalar>
382 bool SolutionState<Scalar>::operator>= (const Scalar& t) const
383 {
384  return (this->metaData_->getTime() >= t);
385 }
386 
387 template<class Scalar>
389 {
390  return (this->metaData_->getTime() == ss.metaData_->getTime());
391 }
392 
393 template<class Scalar>
394 bool SolutionState<Scalar>::operator== (const Scalar& t) const
395 {
396  return (this->metaData_->getTime() == t);
397 }
398 
399 template<class Scalar>
401 {
402  std::string name = "Tempus::SolutionState";
403  return (name);
404 }
405 
406 template<class Scalar>
409  const Teuchos::EVerbosityLevel verbLevel) const
410 {
411  if (verbLevel == Teuchos::VERB_MEDIUM) {
412  out << "(index =" <<std::setw(6)<< this->getIndex()
413  << "; time =" <<std::setw(10)<<std::setprecision(3)<<this->getTime()
414  << "; dt =" <<std::setw(10)<<std::setprecision(3)<<this->getTimeStep()
415  << ")" << std::endl;
416  }
417 
418  if (verbLevel == Teuchos::VERB_EXTREME) {
419  out << description() << "::describe:" << std::endl
420  << "metaData = " << std::endl;
421  metaData_->describe(out,verbLevel);
422  out << "x = " << std::endl;
423  x_->describe(out,verbLevel);
424  if (xdot_ != Teuchos::null) {
425  out << "xdot_ = " << std::endl;
426  xdot_->describe(out,verbLevel);
427  }
428  if (xdotdot_ != Teuchos::null) {
429  out << "xdotdot = " << std::endl;
430  xdotdot_->describe(out,verbLevel);
431  }
432  if (stepperState_ != Teuchos::null) {
433  out << "stepperState = " << std::endl;
434  stepperState_->describe(out,verbLevel);
435  }
436  if (physicsState_ != Teuchos::null) {
437  out << "physicsState = " << std::endl;
438  physicsState_->describe(out,verbLevel);
439  }
440  }
441 }
442 
443 
444 template<class Scalar>
446  const Teuchos::RCP<const SolutionState<Scalar> >& ssIn)
447 {
448  if (!getComputeNorms()) return;
449 
450  auto x = this->getX();
451  this->setXNormL2(Thyra::norm(*x));
452 
453  if (ssIn != Teuchos::null) {
454  auto xIn = ssIn->getX();
455 
456  // dx = x - xIn
457  Teuchos::RCP<Thyra::VectorBase<Scalar> > dx = Thyra::createMember(x->space());
458  Thyra::V_VmV(dx.ptr(), *x, *xIn);
459  Scalar dxNorm = Thyra::norm(*dx);
460  Scalar xInNorm = Thyra::norm(*xIn);
461  this->setDxNormL2Abs(dxNorm);
462  // Compute change, e.g., ||x^n-x^(n-1)||/||x^(n-1)||
463  const Scalar eps = std::numeric_limits<Scalar>::epsilon();
464  const Scalar min = std::numeric_limits<Scalar>::min();
465  if ( xInNorm < min/eps ) { // numerically zero
466  this->setDxNormL2Rel(std::numeric_limits<Scalar>::infinity());
467  } else {
468  //this->setDxNormL2Rel(dxNorm/(xInNorm + eps));
469  this->setDxNormL2Rel(dxNorm/(xInNorm*(1.0 + 1.0e4*eps)));
470  }
471  }
472 }
473 
474 
475 // Nonmember constructors.
476 // ------------------------------------------------------------------------
477 
478 template<class Scalar>
482  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot)
483 {
485  metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
486 
488  stepperState_nc = Teuchos::rcp(new StepperState<Scalar>("Default"));
489 
491  physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ());
492 
494  metaData_nc, x, xdot, xdotdot, stepperState_nc, physicsState_nc));
495 
496  return ss;
497 }
498 
499 template<class Scalar>
501  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
502  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
503  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot)
504 {
507 
509  stepperState = Teuchos::rcp(new StepperState<Scalar>("Default"));
510 
512  physicsState = Teuchos::rcp(new PhysicsState<Scalar> ());
513 
515  metaData, x, xdot, xdotdot, stepperState, physicsState));
516 
517  return ss;
518 }
519 
520 
521 template<class Scalar>
523  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
524  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
525  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
526 {
527  typedef Thyra::ModelEvaluatorBase MEB;
528  using Teuchos::rcp_const_cast;
529 
530  auto metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
531  metaData_nc->setSolutionStatus(Status::PASSED);
532 
533  MEB::InArgs<Scalar> inArgs = model->getNominalValues();
534 
535  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.supports(MEB::IN_ARG_x) == false,
536  std::logic_error,
537  model->description() << "does not support an x solution vector!");
538 
539  // The solution vector, x, is required (usually).
540  auto x_nc = rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x());
541 
542  // The solution derivative, xdot, can be optional provided, based on
543  // application needs. Here we will base it on "supports" IN_ARG_x_dot.
544  // Depending on the stepper used, a temporary xdot vector may be created
545  // within the Stepper, but not moved to the SolutionState.
547  if (inArgs.supports(MEB::IN_ARG_x_dot)) {
548  xdot_nc = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
549  } else {
550  xdot_nc = Teuchos::null;
551  }
552 
553  // Similar as xdot.
555  if (inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
556  xdotdot_nc =
557  rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x_dot_dot());
558  } else {
559  xdotdot_nc = Teuchos::null;
560  }
561 
562  Teuchos::RCP<StepperState<Scalar> > stepperState_nc;
563  if (stepperState == Teuchos::null) {
564  stepperState_nc = Teuchos::rcp(new StepperState<Scalar> ()); // Use default
565  } else {
566  stepperState_nc = stepperState;
567  }
568 
569  Teuchos::RCP<PhysicsState<Scalar> > physicsState_nc;
570  if (physicsState == Teuchos::null) {
571  physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ()); // Use default
572  } else {
573  physicsState_nc = physicsState;
574  }
575 
577  metaData_nc, x_nc, xdot_nc, xdotdot_nc, stepperState_nc, physicsState_nc));
578 
579  return ss;
580 }
581 
582 } // namespace Tempus
583 #endif // Tempus_SolutionState_impl_hpp
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateX(const Teuchos::RCP< Thyra::VectorBase< Scalar > > &x, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdot=Teuchos::null, const Teuchos::RCP< Thyra::VectorBase< Scalar > > &xdotdot=Teuchos::null)
Nonmember constructor from non-const solution vectors, x.
PhysicsState is a simple class to hold information about the physics.
Teuchos::RCP< const SolutionStateMetaData< Scalar > > metaData_
Meta Data for the solution state.
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::RCP< PhysicsState< Scalar > > physicsState_nc_
Teuchos::RCP< const PhysicsState< Scalar > > physicsState_
PhysicsState for this SolutionState.
StepperState is a simple class to hold state information about the stepper.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
bool operator>(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
Teuchos::RCP< const StepperState< Scalar > > stepperState_
StepperState for this SolutionState.
Ptr< T > ptr() const
virtual void copySolutionData(const Teuchos::RCP< const SolutionState< Scalar > > &s)
Deep copy solution data, but keep metaData untouched.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
SolutionState()
Default Constructor – Not meant for immediate adding to SolutionHistory. This constructor does not se...
bool operator>=(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
Teuchos::RCP< SolutionState< Scalar > > createSolutionStateME(const Teuchos::RCP< const Thyra::ModelEvaluator< Scalar > > &model, const Teuchos::RCP< StepperState< Scalar > > &stepperState=Teuchos::null, const Teuchos::RCP< PhysicsState< Scalar > > &physicsState=Teuchos::null)
Nonmember constructor from Thyra ModelEvaluator.
bool operator==(const SolutionState< Scalar > &ss) const
Equality comparison for matching.
Teuchos::RCP< StepperState< Scalar > > stepperState_nc_
bool operator<=(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
virtual Teuchos::RCP< SolutionState< Scalar > > clone() const
This is a deep copy constructor.
virtual void computeNorms(const Teuchos::RCP< const SolutionState< Scalar > > &ssIn=Teuchos::null)
Compute the solution norms, and solution change from ssIn, if provided.
bool operator<(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
virtual std::string description() const
Solution state for integrators and steppers. SolutionState contains the metadata for solutions and th...
Teuchos::RCP< SolutionStateMetaData< Scalar > > metaData_nc_
virtual void copy(const Teuchos::RCP< const SolutionState< Scalar > > &ss)
This is a deep copy.