Skip to content

Commit

Permalink
Tempus: Fix Failing Tests Due to FPE
Browse files Browse the repository at this point in the history
Investigated failing tests from Tempus_*_MPI_1 failing in Trilinos-atdm-ats1-knl_intel-19.0.4_mpich-7.7.15_openmp_static_opt build starting 2021-10-15 #9881

 * Excluded the Intel compiler for FPE testing.
 * Fixed an FPE divide-by-zero in Stepper_ErrorNorm
   - Added unit test to cover it
 * Added some doxygen on Stepper_ErrorNorm
  • Loading branch information
ccober6 committed Nov 13, 2021
1 parent 37e9418 commit 39ce549
Show file tree
Hide file tree
Showing 5 changed files with 154 additions and 27 deletions.
70 changes: 51 additions & 19 deletions packages/tempus/src/Tempus_Stepper_ErrorNorm_decl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -16,37 +16,69 @@
#include "Thyra_VectorSpaceFactoryBase.hpp"
namespace Tempus {

/** \brief Stepper_ErrorNorm provides error norm calcualtions for variable time stepping.
*
*/
template<class Scalar>
class Stepper_ErrorNorm
{
public:

public:
/// Default Constructor
Stepper_ErrorNorm();

// ctor
Stepper_ErrorNorm();
/// Constructor
Stepper_ErrorNorm(const Scalar relTol, const Scalar absTol);

Stepper_ErrorNorm(const Scalar relTol, const Scalar absTol);
/// Destructor
~Stepper_ErrorNorm() {};

~Stepper_ErrorNorm() {};
/** \brief Compute the weigthed root mean square norm.
*
* The WRMS norm is
* \f[
* e_{wrms} \equiv \sqrt{ \frac{1}{N} \sum_{i=1}^N \left( \frac
* {u^n_i}{A_{tol} + \max (|u^n_i|, |u^{n+1}_i|) R_{tol}} \right) ^2 }
* \f]
* where
* - \f$A_{tol}\f$ is the absolute tolerance,
* - \f$R_{tol}\f$ is the relative tolerance,
* - \f$\max\f$ is the pairwise maximum,
* - \f$u^n\f$ is the current solution, and
* - \f$u^{n+1}\f$ is the next solution.
* - \f$ u_i^n\f$ denotes component \f$i\f$ of time step \f$n\f$.
* - \f$ N\f$ denotes the number of unknowns
*/
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);

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);
/** \brief Compute the error Norm.
*
* The error norm is
* \f[
* e = \max_i ( u_i / ( A_{tol} + |u_i| R_{tol}))
* \f]
* where
* \f$A_{tol}\f$ is the absolute tolerance,
* \f$R_{tol}\f$ is the relative tolerance, and
* \f$\max_i\f$ is the maximum over all elements in \f$x\f$.
*/
Scalar errorNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x);

Scalar errorNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x);
void setRelativeTolerance(const Scalar relTol) { relTol_ = relTol; }
void setAbsoluteTolerance(const Scalar absTol) { abssTol_ = absTol; }

void setRelativeTolerance(const Scalar relTol) { relTol_ = relTol; }
void setAbsoluteTolerance(const Scalar absTol) { abssTol_ = absTol; }


protected:
protected:

Scalar relTol_;
Scalar abssTol_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> u_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> uNext_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> errorWeightVector_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> scratchVector_;
Scalar relTol_;
Scalar abssTol_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> u_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> uNext_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> errorWeightVector_;
Teuchos::RCP<Thyra::VectorBase<Scalar>> scratchVector_;

};

Expand Down
26 changes: 21 additions & 5 deletions packages/tempus/src/Tempus_Stepper_ErrorNorm_impl.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,10 @@
#include "Thyra_MultiVectorStdOps_decl.hpp"
#include "Thyra_VectorSpaceBase_decl.hpp"
#include "Thyra_VectorStdOps_decl.hpp"

#include "Tempus_NumericalUtils.hpp"


namespace Tempus {

template<class Scalar>
Expand All @@ -27,7 +31,7 @@ Stepper_ErrorNorm<Scalar>::Stepper_ErrorNorm(const Scalar relTol, const Scalar a

template<class Scalar>
Scalar Stepper_ErrorNorm<Scalar>::
computeWRMSNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x,
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)
{
Expand All @@ -44,7 +48,13 @@ computeWRMSNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x,
Thyra::abs(*x, u_.ptr());
Thyra::abs(*xNext, uNext_.ptr());
Thyra::pair_wise_max_update(relTol_, *u_, uNext_.ptr());
Thyra::add_scalar(abssTol_, 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());
Expand All @@ -59,17 +69,23 @@ template<class Scalar>
Scalar Stepper_ErrorNorm<Scalar>::
errorNorm(const Teuchos::RCP<const Thyra::VectorBase<Scalar>> &x)
{
if (scratchVector_ == Teuchos::null)
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_);
Thyra::Vp_S(scratchVector_.ptr(), abssTol_);
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;

}


Expand Down
6 changes: 3 additions & 3 deletions packages/tempus/test/TestUtils/Tempus_UnitTestMainUtils.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
#ifndef TEMPUS_UNIT_TEST_MAIN_UTILS_HPP
#define TEMPUS_UNIT_TEST_MAIN_UTILS_HPP

#if defined(__linux__) && defined(__GNUC__)
#if defined(__linux__) && defined(__GNUC__) && !defined(__INTEL_COMPILER)
#include <fenv.h>
#elif defined(__APPLE__) && defined(__GNUC__) && defined(__SSE__)
#include <xmmintrin.h>
Expand All @@ -33,7 +33,7 @@ void enableFPE(bool enableFPE)
#endif

if (enableFPE) {
#if defined(__linux__) && defined(__GNUC__)
#if defined(__linux__) && defined(__GNUC__) && !defined(__INTEL_COMPILER)
feenableexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_INVALID);
#elif defined(__APPLE__) && defined(__GNUC__) && defined(__SSE__)
eMask = _MM_GET_EXCEPTION_MASK(); // Save current eMask so we can disable.
Expand All @@ -42,7 +42,7 @@ void enableFPE(bool enableFPE)
& ~_MM_MASK_INVALID);
#endif
} else {
#if defined(__linux__) && defined(__GNUC__)
#if defined(__linux__) && defined(__GNUC__) && !defined(__INTEL_COMPILER)
fedisableexcept(FE_DIVBYZERO | FE_OVERFLOW | FE_INVALID);
#elif defined(__APPLE__) && defined(__GNUC__) && defined(__SSE__)
_MM_SET_EXCEPTION_MASK(eMask);
Expand Down
7 changes: 7 additions & 0 deletions packages/tempus/unit_test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -421,6 +421,13 @@ TRIBITS_ADD_EXECUTABLE_AND_TEST(
NUM_MPI_PROCS 1
)

TRIBITS_ADD_EXECUTABLE_AND_TEST(
UnitTest_Stepper_ErrorNorm
SOURCES Tempus_UnitTest_Stepper_ErrorNorm.cpp ${TEMPUS_UNIT_TEST_MAIN}
TESTONLYLIBS tempus_test_models
NUM_MPI_PROCS 1
)

TRIBITS_COPY_FILES_TO_BINARY_DIR(UnitTest_NewmarkImplicitAForm_CopyFiles
DEST_FILES Tempus_NewmarkImplicitAForm_HarmonicOscillator_Damped_SecondOrder.xml Tempus_NewmarkExplicitAForm_HarmonicOscillator_Damped.xml Tempus_DIRK_VanDerPol.xml
EXEDEPS UnitTest_NewmarkImplicitAForm
Expand Down
72 changes: 72 additions & 0 deletions packages/tempus/unit_test/Tempus_UnitTest_Stepper_ErrorNorm.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
// @HEADER
// ****************************************************************************
// Tempus: Copyright (2017) Sandia Corporation
//
// Distributed under BSD 3-clause license (See accompanying file Copyright.txt)
// ****************************************************************************
// @HEADER

#include "Tempus_UnitTest_Utils.hpp"

#include "Thyra_DefaultSpmdVectorSpace.hpp"

#include "Tempus_NumericalUtils.hpp"
#include "Tempus_Stepper_ErrorNorm.hpp"

#include "../TestModels/DahlquistTestModel.hpp"


namespace Tempus_Unit_Test {

// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(Stepper_ErrorNorm, computeWRMSNorm)
{

int N = 3;
Teuchos::RCP<const Thyra::VectorSpaceBase<double> > xSpace =
Thyra::defaultSpmdVectorSpace<double>(N);
auto x = Thyra::createMember(xSpace);
auto xNext = Thyra::createMember(xSpace);
auto eVec = Thyra::createMember(xSpace);

auto tol = Tempus::numericalTol<double>();
auto eNorm = Teuchos::rcp(new Tempus::Stepper_ErrorNorm<double>(tol, 0.0));

Thyra::assign(x.ptr(), 0.0);
Thyra::assign(xNext.ptr(), 0.0);
Thyra::assign(eVec.ptr(), 0.0);
double rmsNorm = eNorm->computeWRMSNorm(x, xNext, eVec);
TEST_FLOATING_EQUALITY(rmsNorm, 0.0, 1.0e-12);

Thyra::assign(x.ptr(), 1.0);
Thyra::assign(xNext.ptr(), 1.0+tol);
Thyra::assign(eVec.ptr(), tol);
rmsNorm = eNorm->computeWRMSNorm(x, xNext, eVec);
TEST_FLOATING_EQUALITY(rmsNorm, 1.0/std::sqrt(N), 1.0e-12);
}


// ************************************************************
// ************************************************************
TEUCHOS_UNIT_TEST(Stepper_ErrorNorm, errorNorm)
{

Teuchos::RCP<const Thyra::VectorSpaceBase<double> > xSpace =
Thyra::defaultSpmdVectorSpace<double>(3);
auto x = Thyra::createMember(xSpace);

auto tol = Tempus::numericalTol<double>();
auto eNorm = Teuchos::rcp(new Tempus::Stepper_ErrorNorm<double>(tol, tol));

Thyra::assign(x.ptr(), 0.0);
double norm = eNorm->errorNorm(x);
TEST_FLOATING_EQUALITY(norm, tol, 1.0e-12);

Thyra::assign(x.ptr(), 1.0);
norm = eNorm->errorNorm(x);
TEST_FLOATING_EQUALITY(norm, 1.0/(2.0*tol), 1.0e-12);
}


} // namespace Tempus_Test

0 comments on commit 39ce549

Please sign in to comment.