Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
14namespace Tempus {
15
16
17template<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{
30 metaData_nc_ = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
32 stepperState_nc_ = Teuchos::rcp(new StepperState<Scalar>("Default"));
34 physicsState_nc_ = Teuchos::rcp(new PhysicsState<Scalar> ());
36}
37
38
39template<class Scalar>
41 const Teuchos::RCP<SolutionStateMetaData<Scalar> > metaData,
42 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& x,
43 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdot,
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) {
61 stepperState_nc_ = Teuchos::rcp(new StepperState<Scalar>("Default"));
63 }
64 if (physicsState_nc_ == Teuchos::null) {
65 physicsState_nc_ = Teuchos::rcp(new PhysicsState<Scalar> ());
67 }
68}
69
70template<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) {
92 stepperState_ = Teuchos::rcp(new StepperState<Scalar>("Default"));
93 }
94 if (physicsState_ == Teuchos::null) {
95 physicsState_ = Teuchos::rcp(new PhysicsState<Scalar> ());
96 }
97}
98
99
100template<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
117template<class Scalar>
118Teuchos::RCP<SolutionState<Scalar> > SolutionState<Scalar>::clone() const
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
147template<class Scalar>
149copy(const Teuchos::RCP<const SolutionState<Scalar> >& ss)
150{
151 metaData_nc_->copy(ss->metaData_);
152 this->copySolutionData(ss);
153}
154
155
156template<class Scalar>
158copySolutionData(const Teuchos::RCP<const SolutionState<Scalar> >& ss)
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
212template<class Scalar>
214{
215 return (this->metaData_->getTime() < ss.metaData_->getTime());
216}
217
218template<class Scalar>
220{
221 return (this->metaData_->getTime() <= ss.metaData_->getTime());
222}
223
224template<class Scalar>
225bool SolutionState<Scalar>::operator< (const Scalar& t) const
226{
227 return (this->metaData_->getTime() < t);
228}
229
230template<class Scalar>
231bool SolutionState<Scalar>::operator<= (const Scalar& t) const
232{
233 return (this->metaData_->getTime() <= t);
234}
235
236template<class Scalar>
238{
239 return (this->metaData_->getTime() > ss.metaData_->getTime());
240}
241
242template<class Scalar>
244{
245 return (this->metaData_->getTime() >= ss.metaData_->getTime());
246}
247
248template<class Scalar>
249bool SolutionState<Scalar>::operator> (const Scalar& t) const
250{
251 return (this->metaData_->getTime() > t);
252}
253
254template<class Scalar>
255bool SolutionState<Scalar>::operator>= (const Scalar& t) const
256{
257 return (this->metaData_->getTime() >= t);
258}
259
260template<class Scalar>
262{
263 return (this->metaData_->getTime() == ss.metaData_->getTime());
264}
265
266template<class Scalar>
267bool SolutionState<Scalar>::operator== (const Scalar& t) const
268{
269 return (this->metaData_->getTime() == t);
270}
271
272template<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
284template<class Scalar>
286 Teuchos::FancyOStream &out,
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
320template<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
354template<class Scalar>
355Teuchos::RCP<SolutionState<Scalar> > createSolutionStateX(
356 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& x,
357 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdot,
358 const Teuchos::RCP<Thyra::VectorBase<Scalar> >& xdotdot)
359{
360 Teuchos::RCP<SolutionStateMetaData<Scalar> >
361 metaData_nc = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
362
363 Teuchos::RCP<StepperState<Scalar> >
364 stepperState_nc = Teuchos::rcp(new StepperState<Scalar>("Default"));
365
366 Teuchos::RCP<PhysicsState<Scalar> >
367 physicsState_nc = Teuchos::rcp(new PhysicsState<Scalar> ());
368
369 Teuchos::RCP<SolutionState<Scalar> > ss = rcp(new SolutionState<Scalar> (
370 metaData_nc, x, xdot, xdotdot, stepperState_nc, physicsState_nc));
371
372 return ss;
373}
374
375template<class Scalar>
376Teuchos::RCP<SolutionState<Scalar> > createSolutionStateX(
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{
381 Teuchos::RCP<const SolutionStateMetaData<Scalar> >
382 metaData = Teuchos::rcp(new SolutionStateMetaData<Scalar>());
383
384 Teuchos::RCP<const StepperState<Scalar> >
385 stepperState = Teuchos::rcp(new StepperState<Scalar>("Default"));
386
387 Teuchos::RCP<const PhysicsState<Scalar> >
388 physicsState = Teuchos::rcp(new PhysicsState<Scalar> ());
389
390 Teuchos::RCP<SolutionState<Scalar> > ss = rcp(new SolutionState<Scalar> (
391 metaData, x, xdot, xdotdot, stepperState, physicsState));
392
393 return ss;
394}
395
396
397template<class Scalar>
398Teuchos::RCP<SolutionState<Scalar> > createSolutionStateME(
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.
422 Teuchos::RCP<Thyra::VectorBase<Scalar> > xdot_nc;
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.
430 Teuchos::RCP<Thyra::VectorBase<Scalar> > xdotdot_nc;
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
452 Teuchos::RCP<SolutionState<Scalar> > ss = rcp(new SolutionState<Scalar> (
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
PhysicsState is a simple class to hold information about the physics.
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.
virtual void computeNorms(const Teuchos::RCP< const SolutionState< Scalar > > &ssIn=Teuchos::null)
Compute the solution norms, and solution change from ssIn, if provided.
Teuchos::RCP< StepperState< Scalar > > stepperState_nc_
Teuchos::RCP< const SolutionStateMetaData< Scalar > > metaData_
Meta Data for the solution state.
virtual std::string description() const
virtual void copySolutionData(const Teuchos::RCP< const SolutionState< Scalar > > &s)
Deep copy solution data, but keep metaData untouched.
SolutionState()
Default Constructor – Not meant for immediate adding to SolutionHistory. This constructor does not se...
Teuchos::RCP< const PhysicsState< Scalar > > physicsState_
PhysicsState for this SolutionState.
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
bool operator==(const SolutionState< Scalar > &ss) const
Equality comparison for matching.
bool operator<(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
bool operator>(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
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.
bool operator>=(const SolutionState< Scalar > &ss) const
Less than comparison for sorting based on time.
Teuchos::RCP< PhysicsState< Scalar > > physicsState_nc_
virtual Teuchos::RCP< SolutionState< Scalar > > clone() const
This is a deep copy constructor.
StepperState is a simple class to hold state information about the stepper.
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.
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.