add FPEs to mam4xx unit testing (#479)

* add FPEs to mam4xx unit testing

* math_tests--address corner case in root_finders.hpp newton solver for when derivative is zero at initial guess

* clang-format

haero/root_finders.hpp

@@ -59,9 +59,18 @@ template <typename ScalarFunction> struct NewtonSolver {
                               value_type>::type
       solve_impl() {
     bool keep_going = true;
+    value_type xnp1;
+    value_type f_deriv = 0;
     while (keep_going) {
       ++counter;
-      const value_type xnp1 = xroot - f(xroot) / f.derivative(xroot);
+      f_deriv = f.derivative(xroot);
+      if (FloatingPoint<value_type>::zero(f_deriv)) {
+        xroot = nan("");
+        keep_going = false;
+        fail = true;
+        break;
+      }
+      xnp1 = xroot - f(xroot) / f_deriv;
       iter_diff = abs(xnp1 - xroot);
       keep_going = !(FloatingPoint<value_type>::zero(iter_diff, conv_tol));
       EKAT_KERNEL_ASSERT_MSG(counter <= max_iter,

haero/testing.cpp

@@ -128,6 +128,9 @@ Surface create_surface() { return Surface(); }
 void ekat_initialize_test_session(int argc, char **argv,
                                   const bool print_config) {
   ekat::initialize_ekat_session(argc, argv, print_config);
+#ifdef EKAT_ENABLE_FPE
+  ekat::enable_fpes(haero::testing::default_fpes);
+#endif
 }
 
 // This implementation of ekat_finalize_test_session calls

haero/testing.hpp

@@ -5,6 +5,8 @@
 #ifndef HAERO_TESTING_HPP
 #define HAERO_TESTING_HPP
 
+#include <cfenv>
+
 #include <haero/atmosphere.hpp>
 #include <haero/surface.hpp>
 
@@ -14,6 +16,8 @@ namespace haero {
 // environments.
 namespace testing {
 
+constexpr int default_fpes = FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW;
+
 /// creates an Atmosphere object that stores a column of data with the given
 /// number of vertical levels and the given planetary boundary height
 /// @param [in] num_levels the number of vertical levels per column stored by

haero/tests/math_tests.hpp

@@ -105,6 +105,7 @@ template <typename T> struct MonicParabola {
   Real b;
 
   KOKKOS_INLINE_FUNCTION
+  // y = x**2 + 3
   MonicParabola() : a(0), b(3) {}
 
   KOKKOS_INLINE_FUNCTION