Tempus  Version of the Day
Time Integration
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Tempus_Stepper_ErrorNorm_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_Stepper_ErrorNorm_impl_hpp
10 #define Tempus_Stepper_ErrorNorm_impl_hpp
11 
13 #include "Thyra_DefaultSerialDenseLinearOpWithSolve_decl.hpp"
14 #include "Thyra_MultiVectorStdOps_decl.hpp"
15 #include "Thyra_VectorSpaceBase_decl.hpp"
16 #include "Thyra_VectorStdOps_decl.hpp"
17 
19 
20 
21 namespace Tempus {
22 
23 template<class Scalar>
24 Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm() : relTol_(1.0e-4), abssTol_(1.0e-4)
25 {}
26 
27 template<class Scalar>
28 Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm(const Scalar relTol, const Scalar absTol) :
29  relTol_(relTol), abssTol_(absTol)
30 {}
31 
32 template<class Scalar>
35  const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &xNext,
36  const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &err)
37 {
38  if (errorWeightVector_ == Teuchos::null)
39  errorWeightVector_ = Thyra::createMember(x->space());
40 
41  if (u_ == Teuchos::null)
42  u_ = Thyra::createMember(x->space());
43 
44  if (uNext_ == Teuchos::null)
45  uNext_ = Thyra::createMember(x->space());
46 
47  // Compute: Atol + max(|u^n|, |u^{n+1}| ) * Rtol
48  Thyra::abs(*x, u_.ptr());
49  Thyra::abs(*xNext, uNext_.ptr());
50  Thyra::pair_wise_max_update(relTol_, *u_, uNext_.ptr());
51  if ( !approxZero(abssTol_) ) {
52  Thyra::add_scalar(abssTol_, uNext_.ptr());
53  } else {
54  Scalar absTol = Thyra::norm_2(*uNext_)*numericalTol<Scalar>();
55  if ( approxZero(absTol) ) absTol = numericalTol<Scalar>();
56  Thyra::add_scalar(absTol, uNext_.ptr());
57  }
58 
59  Thyra::assign(errorWeightVector_.ptr(), Teuchos::ScalarTraits<Scalar>::zero());
60  Thyra::ele_wise_divide(Teuchos::as<Scalar>(1.0), *err, *uNext_, errorWeightVector_.ptr());
61 
62  const auto space_dim = err->space()->dim();
63  Scalar err_norm = std::abs( Thyra::norm(*errorWeightVector_) / space_dim);
64  return err_norm;
65 }
66 
67 
68 template<class Scalar>
71 {
72  if (scratchVector_ == Teuchos::null)
73  scratchVector_ = Thyra::createMember(x->space());
74 
75  Thyra::assign(scratchVector_.ptr(), *x); // | U |
76  Thyra::abs(*x, scratchVector_.ptr());
77  Thyra::Vt_S(scratchVector_.ptr(), relTol_);
78  if ( !approxZero(abssTol_) ) {
79  Thyra::Vp_S(scratchVector_.ptr(), abssTol_);
80  } else {
81  Scalar absTol = Thyra::norm_2(*scratchVector_)*numericalTol<Scalar>();
82  if ( approxZero(absTol) ) absTol = numericalTol<Scalar>();
83  Thyra::add_scalar(absTol, scratchVector_.ptr());
84  }
85  Thyra::ele_wise_divide(Teuchos::as<Scalar>(1.0), *x, *scratchVector_, scratchVector_.ptr());
86  Scalar err = Thyra::norm_inf(*scratchVector_);
87  return err;
88 
89 }
90 
91 
92 }
93 
94 #endif
Scalar computeWRMSNorm(const Teuchos::RCP< const Thyra::VectorBase< Scalar >> &x, const Teuchos::RCP< const Thyra::VectorBase< Scalar >> &xNext, const Teuchos::RCP< const Thyra::VectorBase< Scalar >> &err)
Compute the weigthed root mean square norm.
bool approxZero(Scalar value, Scalar tol=Teuchos::ScalarTraits< Scalar >::sfmin())
Test if value is approximately zero within tolerance.
Scalar errorNorm(const Teuchos::RCP< const Thyra::VectorBase< Scalar >> &x)
Compute the error Norm.