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 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  : metaData_ (metaData),
48  metaData_nc_ (metaData),
49  x_ (x),
50  x_nc_ (x),
51  xdot_ (xdot),
52  xdot_nc_ (xdot),
53  xdotdot_ (xdotdot),
54  xdotdot_nc_ (xdotdot),
55  stepperState_ (stepperState),
56  stepperState_nc_(stepperState),
57  physicsState_ (physicsState),
58  physicsState_nc_(physicsState)
59 {
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 SolutionStateMetaData<Scalar> > metaData,
73  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
74  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
75  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot,
76  const Teuchos::RCP<const StepperState<Scalar> >& stepperState,
77  const Teuchos::RCP<const PhysicsState<Scalar> >& physicsState)
78  : metaData_ (metaData),
79  metaData_nc_ (Teuchos::null),
80  x_ (x),
81  x_nc_ (Teuchos::null),
82  xdot_ (xdot),
83  xdot_nc_ (Teuchos::null),
84  xdotdot_ (xdotdot),
85  xdotdot_nc_ (Teuchos::null),
86  stepperState_ (stepperState),
87  stepperState_nc_(Teuchos::null),
88  physicsState_ (physicsState),
89  physicsState_nc_(Teuchos::null)
90 {
91  if (stepperState_ == Teuchos::null) {
93  }
94  if (physicsState_ == Teuchos::null) {
96  }
97 }
98 
99 
100 template<class Scalar>
102  :metaData_ (ss_.metaData_),
103  metaData_nc_ (ss_.metaData_nc_),
104  x_ (ss_.x_),
105  x_nc_ (ss_.x_nc_),
106  xdot_ (ss_.xdot_),
107  xdot_nc_ (ss_.xdot_nc_),
108  xdotdot_ (ss_.xdotdot_),
109  xdotdot_nc_ (ss_.xdotdot_nc_),
110  stepperState_ (ss_.stepperState_),
111  stepperState_nc_(ss_.stepperState_nc_),
112  physicsState_ (ss_.physicsState_),
113  physicsState_nc_(ss_.physicsState_nc_)
114 {}
115 
116 
117 template<class Scalar>
119 {
120  using Teuchos::RCP;
121 
122  RCP<SolutionStateMetaData<Scalar> > metaData_out;
123  if (!Teuchos::is_null(metaData_)) metaData_out = metaData_->clone();
124 
125  RCP<Thyra::VectorBase<Scalar> > x_out;
126  if (!Teuchos::is_null(x_)) x_out = x_->clone_v();
127 
128  RCP<Thyra::VectorBase<Scalar> > xdot_out;
129  if (!Teuchos::is_null(xdot_)) xdot_out = xdot_->clone_v();
130 
131  RCP<Thyra::VectorBase<Scalar> > xdotdot_out;
132  if (!Teuchos::is_null(xdotdot_)) xdotdot_out = xdotdot_->clone_v();
133 
134  RCP<StepperState<Scalar> > sS_out;
135  if (!Teuchos::is_null(stepperState_)) sS_out=stepperState_->clone();
136 
137  RCP<PhysicsState<Scalar> > pS_out;
138  if (!Teuchos::is_null(physicsState_)) pS_out=physicsState_->clone();
139 
140  RCP<SolutionState<Scalar> > ss_out = Teuchos::rcp(new SolutionState<Scalar> (
141  metaData_out, x_out, xdot_out, xdotdot_out, sS_out, pS_out));
142 
143  return ss_out;
144 }
145 
146 
147 template<class Scalar>
150 {
151  metaData_nc_->copy(ss->metaData_);
152  this->copySolutionData(ss);
153 }
154 
155 
156 template<class Scalar>
159 {
160  if (ss->x_ == Teuchos::null)
161  x_nc_ = Teuchos::null;
162  else {
163  if (x_nc_ == Teuchos::null) {
164  x_nc_ = ss->x_->clone_v();
165  }
166  else
167  Thyra::V_V(x_nc_.ptr(), *(ss->x_));
168  }
169  x_ = x_nc_;
170 
171  if (ss->xdot_ == Teuchos::null)
172  xdot_nc_ = Teuchos::null;
173  else {
174  if (xdot_nc_ == Teuchos::null)
175  xdot_nc_ = ss->xdot_->clone_v();
176  else
177  Thyra::V_V(xdot_nc_.ptr(), *(ss->xdot_));
178  }
179  xdot_ = xdot_nc_;
180 
181  if (ss->xdotdot_ == Teuchos::null)
182  xdotdot_nc_ = Teuchos::null;
183  else {
184  if (xdotdot_nc_ == Teuchos::null)
185  xdotdot_nc_ = ss->xdotdot_->clone_v();
186  else
187  Thyra::V_V(xdotdot_nc_.ptr(), *(ss->xdotdot_));
188  }
189  xdotdot_ = xdotdot_nc_;
190 
191  if (ss->stepperState_ == Teuchos::null)
192  stepperState_nc_ = Teuchos::null;
193  else {
194  if (stepperState_nc_ == Teuchos::null)
195  stepperState_nc_ = ss->stepperState_->clone();
196  else
197  stepperState_nc_->copy(ss->stepperState_);
198  }
199  stepperState_ = stepperState_nc_;
200 
201  if (ss->physicsState_ == Teuchos::null)
202  physicsState_nc_ = Teuchos::null;
203  else {
204  if (physicsState_nc_ == Teuchos::null)
205  physicsState_nc_ = ss->physicsState_->clone();
206  else
207  physicsState_nc_->copy(ss->physicsState_);
208  }
209  physicsState_ = physicsState_nc_;
210 }
211 
212 template<class Scalar>
214 {
215  return (this->metaData_->getTime() < ss.metaData_->getTime());
216 }
217 
218 template<class Scalar>
220 {
221  return (this->metaData_->getTime() <= ss.metaData_->getTime());
222 }
223 
224 template<class Scalar>
225 bool SolutionState<Scalar>::operator< (const Scalar& t) const
226 {
227  return (this->metaData_->getTime() < t);
228 }
229 
230 template<class Scalar>
231 bool SolutionState<Scalar>::operator<= (const Scalar& t) const
232 {
233  return (this->metaData_->getTime() <= t);
234 }
235 
236 template<class Scalar>
238 {
239  return (this->metaData_->getTime() > ss.metaData_->getTime());
240 }
241 
242 template<class Scalar>
244 {
245  return (this->metaData_->getTime() >= ss.metaData_->getTime());
246 }
247 
248 template<class Scalar>
249 bool SolutionState<Scalar>::operator> (const Scalar& t) const
250 {
251  return (this->metaData_->getTime() > t);
252 }
253 
254 template<class Scalar>
255 bool SolutionState<Scalar>::operator>= (const Scalar& t) const
256 {
257  return (this->metaData_->getTime() >= t);
258 }
259 
260 template<class Scalar>
262 {
263  return (this->metaData_->getTime() == ss.metaData_->getTime());
264 }
265 
266 template<class Scalar>
267 bool SolutionState<Scalar>::operator== (const Scalar& t) const
268 {
269  return (this->metaData_->getTime() == t);
270 }
271 
272 template<class Scalar>
274 {
275  std::ostringstream out;
276  out << "SolutionState"
277  << " (index =" <<std::setw(6)<< this->getIndex()
278  << "; time =" <<std::setw(10)<<std::setprecision(3)<<this->getTime()
279  << "; dt =" <<std::setw(10)<<std::setprecision(3)<<this->getTimeStep()
280  << ")";
281  return out.str();
282 }
283 
284 template<class Scalar>
287  const Teuchos::EVerbosityLevel verbLevel) const
288 {
289  auto l_out = Teuchos::fancyOStream( out.getOStream() );
290  Teuchos::OSTab ostab(*l_out, 2, this->description());
291  l_out->setOutputToRootOnly(0);
292 
293  *l_out << "\n--- " << this->description() << " ---" << std::endl;
294 
295  if (Teuchos::as<int>(verbLevel) >= Teuchos::as<int>(Teuchos::VERB_EXTREME)) {
296 
297  metaData_->describe(*l_out,verbLevel);
298  *l_out << " x = " << std::endl;
299  x_->describe(*l_out,verbLevel);
300 
301  if (xdot_ != Teuchos::null) {
302  *l_out << " xdot_ = " << std::endl;
303  xdot_->describe(*l_out,verbLevel);
304  }
305  if (xdotdot_ != Teuchos::null) {
306  *l_out << " xdotdot = " << std::endl;
307  xdotdot_->describe(*l_out,verbLevel);
308  }
309 
310  if (stepperState_ != Teuchos::null)
311  stepperState_->describe(*l_out,verbLevel);
312  if (physicsState_ != Teuchos::null)
313  physicsState_->describe(*l_out,verbLevel);
314 
315  *l_out << std::string(this->description().length()+8, '-') <<std::endl;
316  }
317 }
318 
319 
320 template<class Scalar>
322  const Teuchos::RCP<const SolutionState<Scalar> >& ssIn)
323 {
324  if (!getComputeNorms()) return;
325 
326  auto x = this->getX();
327  this->setXNormL2(Thyra::norm(*x));
328 
329  if (ssIn != Teuchos::null) {
330  auto xIn = ssIn->getX();
331 
332  // dx = x - xIn
333  Teuchos::RCP<Thyra::VectorBase<Scalar> > dx = Thyra::createMember(x->space());
334  Thyra::V_VmV(dx.ptr(), *x, *xIn);
335  Scalar dxNorm = Thyra::norm(*dx);
336  Scalar xInNorm = Thyra::norm(*xIn);
337  this->setDxNormL2Abs(dxNorm);
338  // Compute change, e.g., ||x^n-x^(n-1)||/||x^(n-1)||
339  const Scalar eps = std::numeric_limits<Scalar>::epsilon();
340  const Scalar min = std::numeric_limits<Scalar>::min();
341  if ( xInNorm < min/eps ) { // numerically zero
342  this->setDxNormL2Rel(std::numeric_limits<Scalar>::infinity());
343  } else {
344  //this->setDxNormL2Rel(dxNorm/(xInNorm + eps));
345  this->setDxNormL2Rel(dxNorm/(xInNorm*(1.0 + 1.0e4*eps)));
346  }
347  }
348 }
349 
350 
351 // Nonmember constructors.
352 // ------------------------------------------------------------------------
353 
354 template<class Scalar>
358  const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot)
359 {
361  metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
362 
364  stepperState_nc = Teuchos::rcp(new StepperState<Scalar>("Default"));
365 
367  physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ());
368 
370  metaData_nc, x, xdot, xdotdot, stepperState_nc, physicsState_nc));
371 
372  return ss;
373 }
374 
375 template<class Scalar>
377  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& x,
378  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdot,
379  const Teuchos::RCP<const Thyra::VectorBase<Scalar> >& xdotdot)
380 {
383 
385  stepperState = Teuchos::rcp(new StepperState<Scalar>("Default"));
386 
388  physicsState = Teuchos::rcp(new PhysicsState<Scalar> ());
389 
391  metaData, x, xdot, xdotdot, stepperState, physicsState));
392 
393  return ss;
394 }
395 
396 
397 template<class Scalar>
399  const Teuchos::RCP<const Thyra::ModelEvaluator<Scalar> >& model,
400  const Teuchos::RCP<StepperState<Scalar> >& stepperState,
401  const Teuchos::RCP<PhysicsState<Scalar> >& physicsState)
402 {
403  typedef Thyra::ModelEvaluatorBase MEB;
404  using Teuchos::rcp_const_cast;
405 
406  auto metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
407  metaData_nc->setSolutionStatus(Status::PASSED);
408 
409  MEB::InArgs<Scalar> inArgs = model->getNominalValues();
410 
411  TEUCHOS_TEST_FOR_EXCEPTION(inArgs.supports(MEB::IN_ARG_x) == false,
412  std::logic_error,
413  model->description() << "does not support an x solution vector!");
414 
415  // The solution vector, x, is required (usually).
416  auto x_nc = rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x());
417 
418  // The solution derivative, xdot, can be optional provided, based on
419  // application needs. Here we will base it on "supports" IN_ARG_x_dot.
420  // Depending on the stepper used, a temporary xdot vector may be created
421  // within the Stepper, but not moved to the SolutionState.
423  if (inArgs.supports(MEB::IN_ARG_x_dot)) {
424  xdot_nc = rcp_const_cast<Thyra::VectorBase<Scalar> >(inArgs.get_x_dot());
425  } else {
426  xdot_nc = Teuchos::null;
427  }
428 
429  // Similar as xdot.
431  if (inArgs.supports(MEB::IN_ARG_x_dot_dot)) {
432  xdotdot_nc =
433  rcp_const_cast<Thyra::VectorBase<Scalar> > (inArgs.get_x_dot_dot());
434  } else {
435  xdotdot_nc = Teuchos::null;
436  }
437 
438  Teuchos::RCP<StepperState<Scalar> > stepperState_nc;
439  if (stepperState == Teuchos::null) {
440  stepperState_nc = Teuchos::rcp(new StepperState<Scalar> ()); // Use default
441  } else {
442  stepperState_nc = stepperState;
443  }
444 
445  Teuchos::RCP<PhysicsState<Scalar> > physicsState_nc;
446  if (physicsState == Teuchos::null) {
447  physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ()); // Use default
448  } else {
449  physicsState_nc = physicsState;
450  }
451 
453  metaData_nc, x_nc, xdot_nc, xdotdot_nc, stepperState_nc, physicsState_nc));
454 
455  return ss;
456 }
457 
458 } // namespace Tempus
459 #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
Greater 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
Greater than or equal to 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_
RCP< std::basic_ostream< char_type, traits_type > > getOStream()
bool operator<=(const SolutionState< Scalar > &ss) const
Less than or equal to 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.
Teuchos::RCP< SolutionStateMetaData< Scalar > > metaData_nc_
virtual void copy(const Teuchos::RCP< const SolutionState< Scalar > > &ss)
This is a deep copy.