Tempus Version of the Day
Time Integration
Loading...
Searching...
No Matches
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
12#include "Teuchos_ScalarTraitsDecl.hpp"
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
21namespace Tempus {
22
23template<class Scalar>
24Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm() : relTol_(1.0e-4), abssTol_(1.0e-4)
25{}
26
27template<class Scalar>
28Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm(const Scalar relTol, const Scalar absTol) :
29 relTol_(relTol), abssTol_(absTol)
30{}
31
32template<class Scalar>
34computeWRMSNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x,
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
68template<class Scalar>
70errorNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x)
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.
Scalar errorNorm(const Teuchos::RCP< const Thyra::VectorBase< Scalar > > &x)
Compute the error Norm.
bool approxZero(Scalar value, Scalar tol=Teuchos::ScalarTraits< Scalar >::sfmin())
Test if value is approximately zero within tolerance.