Tempus: Fix Failing Tests Due to FPE

packages/tempus/src/Tempus_Stepper_ErrorNorm_decl.hpp

@@ -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_;
 
 };
 

packages/tempus/src/Tempus_Stepper_ErrorNorm_impl.hpp

@@ -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>
@@ -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)
 {
@@ -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());
@@ -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;
-  
+
 }
 
 

packages/tempus/test/TestUtils/Tempus_UnitTestMainUtils.hpp

@@ -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>
@@ -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.
@@ -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);

packages/tempus/unit_test/CMakeLists.txt

@@ -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

packages/tempus/unit_test/Tempus_UnitTest_Stepper_ErrorNorm.cpp

@@ -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