Tempus_Stepper_ErrorNorm_impl.hpp

Go to the documentation of this file.

// @HEADER
// ****************************************************************************
//                Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#ifndef Tempus_Stepper_ErrorNorm_impl_hpp
#define Tempus_Stepper_ErrorNorm_impl_hpp

#include "Teuchos_ScalarTraitsDecl.hpp"
#include "Thyra_DefaultSerialDenseLinearOpWithSolve_decl.hpp"
#include "Thyra_MultiVectorStdOps_decl.hpp"
#include "Thyra_VectorSpaceBase_decl.hpp"
#include "Thyra_VectorStdOps_decl.hpp"

#include "Tempus_NumericalUtils.hpp"

namespace Tempus {

template <class Scalar>
Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm()
  : relTol_(1.0e-4), abssTol_(1.0e-4)
{
}

template <class Scalar>
Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm(const Scalar relTol,
                                             const Scalar absTol)
  : relTol_(relTol), abssTol_(absTol)
{
}

template <class Scalar>
Scalar Stepper_ErrorNorm<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)
{
  if (errorWeightVector_ == Teuchos::null)
    errorWeightVector_ = Thyra::createMember(x->space());

  if (u_ == Teuchos::null) u_ = Thyra::createMember(x->space());

  if (uNext_ == Teuchos::null) uNext_ = Thyra::createMember(x->space());

  // Compute: Atol + max(|u^n|, |u^{n+1}| ) * Rtol
  Thyra::abs(*x, u_.ptr());
  Thyra::abs(*xNext, uNext_.ptr());
  Thyra::pair_wise_max_update(relTol_, *u_, uNext_.ptr());
  if (!approxZero(abssTol_)) {
    Thyra::add_scalar(abssTol_, uNext_.ptr());
  }
  else {
    Scalar absTol = Thyra::norm_2(*uNext_) * numericalTol<Scalar>();
    if (approxZero(absTol)) absTol = numericalTol<Scalar>();
    Thyra::add_scalar(absTol, uNext_.ptr());
  }

  Thyra::assign(errorWeightVector_.ptr(),
                Teuchos::ScalarTraits<Scalar>::zero());
  Thyra::ele_wise_divide(Teuchos::as<Scalar>(1.0), *err, *uNext_,
                         errorWeightVector_.ptr());

  const auto space_dim = err->space()->dim();
  Scalar err_norm      = std::abs(Thyra::norm(*errorWeightVector_) / space_dim);
  return err_norm;
}

template <class Scalar>
Scalar Stepper_ErrorNorm<Scalar>::errorNorm(
    const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x)
{
  if (scratchVector_ == Teuchos::null)
    scratchVector_ = Thyra::createMember(x->space());

  Thyra::assign(scratchVector_.ptr(), *x);  // | U |
  Thyra::abs(*x, scratchVector_.ptr());
  Thyra::Vt_S(scratchVector_.ptr(), relTol_);
  if (!approxZero(abssTol_)) {
    Thyra::Vp_S(scratchVector_.ptr(), abssTol_);
  }
  else {
    Scalar absTol = Thyra::norm_2(*scratchVector_) * numericalTol<Scalar>();
    if (approxZero(absTol)) absTol = numericalTol<Scalar>();
    Thyra::add_scalar(absTol, scratchVector_.ptr());
  }
  Thyra::ele_wise_divide(Teuchos::as<Scalar>(1.0), *x, *scratchVector_,
                         scratchVector_.ptr());
  Scalar err = Thyra::norm_inf(*scratchVector_);
  return err;
}

}  // namespace Tempus

#endif