EAMxx: fix numerical value of fill value

components/eamxx/src/diagnostics/aodvis.cpp

@@ -47,7 +47,7 @@ void AODVis::initialize_impl(const RunType /*run_type*/) {
   auto nondim = ekat::units::Units::nondimensional();
   const auto &grid_name =
       m_diagnostic_output.get_header().get_identifier().get_grid_name();
-  const auto var_fill_value = constants::DefaultFillValue<Real>().value;
+  const auto var_fill_value = constants::fill_value<Real>;
 
   m_mask_val = m_params.get<double>("mask_value", var_fill_value);
 

components/eamxx/src/diagnostics/atm_backtend.cpp

@@ -63,7 +63,7 @@ void AtmBackTendDiag::init_timestep(const util::TimeStamp &start_of_step) {
 }
 
 void AtmBackTendDiag::compute_diagnostic_impl() {
-  Real var_fill_value = constants::DefaultFillValue<Real>().value;
+  Real var_fill_value = constants::fill_value<Real>;
   std::int64_t dt;
 
   const auto &f       = get_field_in(m_name);

components/eamxx/src/diagnostics/field_at_pressure_level.cpp

@@ -85,7 +85,7 @@ initialize_impl (const RunType /*run_type*/)
   // Add a field representing the mask as extra data to the diagnostic field.
   auto nondim = ekat::units::Units::nondimensional();
   const auto& gname = fid.get_grid_name();
-  m_mask_val = m_params.get<double>("mask_value",Real(constants::DefaultFillValue<float>::value));
+  m_mask_val = m_params.get<double>("mask_value",Real(constants::fill_value<double>));
 
 
   std::string mask_name = m_diag_name + " mask";

components/eamxx/src/diagnostics/tests/aodvis_test.cpp

@@ -31,7 +31,7 @@ TEST_CASE("aodvis") {
   using namespace ShortFieldTagsNames;
   using namespace ekat::units;
 
-  Real var_fill_value = constants::DefaultFillValue<Real>().value;
+  Real var_fill_value = constants::fill_value<Real>;
 
   Real some_limit = 0.0025;
 

components/eamxx/src/diagnostics/tests/atm_backtend_test.cpp

@@ -65,7 +65,7 @@ TEST_CASE("atm_backtend") {
   REQUIRE_THROWS(diag_factory.create("AtmBackTendDiag", comm,
                                      params));  // No 'tendency_name'
 
-  Real var_fill_value = constants::DefaultFillValue<Real>().value;
+  Real var_fill_value = constants::fill_value<Real>;
 
   // Set time for qc and randomize its values
   qc.get_header().get_tracking().update_time_stamp(t0);

components/eamxx/src/physics/nudging/eamxx_nudging_process_interface.cpp

@@ -344,7 +344,7 @@ void Nudging::run_impl (const double dt)
     const auto fl = f.get_header().get_identifier().get_layout();
     const auto v  = f.get_view<Real**>();
 
-    Real var_fill_value = constants::DefaultFillValue<Real>().value;
+    Real var_fill_value = constants::fill_value<Real>;
     // Query the helper field for the fill value, if not present use default
     if (f.get_header().has_extra_data("mask_value")) {
       var_fill_value = f.get_header().get_extra_data<Real>("mask_value");

components/eamxx/src/share/io/scorpio_output.hpp

@@ -219,7 +219,7 @@ class AtmosphereOutput
   //   NetCDF: Numeric conversion not representable
   // Also, by default, don't pick max float, to avoid any overflow if the value
   // is used inside other calculation and/or remap.
-  float m_fill_value = constants::DefaultFillValue<float>().value;
+  float m_fill_value = constants::fill_value<float>;
 
   bool m_add_time_dim;
   bool m_track_avg_cnt = false;

components/eamxx/src/share/io/tests/io_basic.cpp

@@ -274,7 +274,7 @@ void read (const std::string& avg_type, const std::string& freq_units,
   // Check that the expected metadata was appropriately set for each variable
   for (const auto& fn: fnames) {
     auto att_fill = scorpio::get_attribute<float>(filename,fn,"_FillValue");
-    REQUIRE(att_fill==constants::DefaultFillValue<float>().value);
+    REQUIRE(att_fill==constants::fill_value<Real>);
 
     auto att_str = scorpio::get_attribute<std::string>(filename,fn,"test");
     REQUIRE (att_str==fn);

components/eamxx/src/share/io/tests/io_filled.cpp

@@ -27,7 +27,7 @@
 namespace scream {
 
 constexpr int num_output_steps = 5;
-constexpr Real FillValue = constants::DefaultFillValue<float>().value;
+constexpr Real FillValue = constants::fill_value<Real>;
 constexpr Real fill_threshold = 0.5;
 
 void set (const Field& f, const double v) {
@@ -255,7 +255,7 @@ void read (const std::string& avg_type, const std::string& freq_units,
   for (const auto& fn: fnames) {
     // NOTE: use float, since default fp_precision for I/O is 'single'
     auto att_fill = scorpio::get_attribute<float>(filename,fn,"_FillValue");
-    REQUIRE(att_fill==constants::DefaultFillValue<float>().value);
+    REQUIRE(att_fill==constants::fill_value<Real>);
   }
 }
 

components/eamxx/src/share/io/tests/io_remap_test.cpp

@@ -283,7 +283,7 @@ TEST_CASE("io_remap_test","io_remap_test")
   {
     // Note, the vertical remapper defaults to a mask value of std numeric limits scaled by 0.1;
     const float mask_val = vert_remap_control.isParameter("Fill Value")
-                         ? vert_remap_control.get<double>("Fill Value") : constants::DefaultFillValue<float>().value;
+                         ? vert_remap_control.get<double>("Fill Value") : constants::fill_value<Real>;
     print ("    -> vertical remap ... \n",io_comm);
     auto gm_vert   = get_test_gm(io_comm,ncols_src,nlevs_tgt);
     auto grid_vert = gm_vert->get_grid("point_grid");
@@ -353,7 +353,7 @@ TEST_CASE("io_remap_test","io_remap_test")
   {
     // Note, the vertical remapper defaults to a mask value of std numeric limits scaled by 0.1;
     const float mask_val = horiz_remap_control.isParameter("Fill Value")
-                         ? horiz_remap_control.get<double>("Fill Value") : constants::DefaultFillValue<float>().value;
+                         ? horiz_remap_control.get<double>("Fill Value") : constants::fill_value<Real>;
     print ("    -> horizontal remap ... \n",io_comm);
     auto gm_horiz   = get_test_gm(io_comm,ncols_tgt,nlevs_src);
     auto grid_horiz = gm_horiz->get_grid("point_grid");
@@ -439,7 +439,7 @@ TEST_CASE("io_remap_test","io_remap_test")
   //                                ---  Vertical + Horizontal Remapping ---
   {
     const float mask_val = vert_horiz_remap_control.isParameter("Fill Value")
-                         ? vert_horiz_remap_control.get<double>("Fill Value") : constants::DefaultFillValue<float>().value;
+                         ? vert_horiz_remap_control.get<double>("Fill Value") : constants::fill_value<Real>;
     print ("    -> vertical + horizontal remap ... \n",io_comm);
     auto gm_vh   = get_test_gm(io_comm,ncols_tgt,nlevs_tgt);
     auto grid_vh = gm_vh->get_grid("point_grid");

components/eamxx/src/share/io/tests/output_restart.cpp

@@ -27,7 +27,7 @@
 
 namespace scream {
 
-constexpr Real FillValue = constants::DefaultFillValue<float>().value;
+constexpr Real FillValue = constants::fill_value<Real>;
 
 std::shared_ptr<FieldManager>
 get_test_fm(const std::shared_ptr<const AbstractGrid>& grid);

components/eamxx/src/share/tests/CMakeLists.txt

@@ -5,6 +5,9 @@ if (NOT SCREAM_ONLY_GENERATE_BASELINES)
   # Test utils
   CreateUnitTest(utils "utils_tests.cpp")
 
+  # Test combine operations
+  CreateUnitTest(combine_ops "combine_ops.cpp")
+
   # Test column ops
   CreateUnitTest(column_ops "column_ops.cpp")
 

components/eamxx/src/share/tests/column_ops.cpp

@@ -7,37 +7,6 @@
 
 namespace {
 
-TEST_CASE ("combine_ops") {
-  using namespace scream;
-  using pack_type = ekat::Pack<Real,SCREAM_PACK_SIZE>;
-
-  constexpr auto Replace      = CombineMode::Replace;
-  constexpr auto Update       = CombineMode::Update;
-  constexpr auto Multiply     = CombineMode::Multiply;
-  constexpr auto Divide       = CombineMode::Divide;
-
-  const pack_type two (2.0);
-  const pack_type four (4.0);
-  const pack_type six (6.0);
-  const pack_type ten (10.0);
-  pack_type x;
-
-  x = two;
-  combine<Replace>(two,x,1,0);
-  REQUIRE ( (x==two).all() );
-
-  combine<Update>(two,x,2.0,1.0);
-  REQUIRE ( (x==six).all() );
-
-  x = two;
-  combine<Multiply>(two,x,1,1);
-  REQUIRE ( (x==four).all() );
-
-  x = four;
-  combine<Divide>(two,x,1,1);
-  REQUIRE ( (x==two).all() );
-}
-
 TEST_CASE("column_ops_ps_1") {
   using namespace scream;
   using device_type = DefaultDevice;

components/eamxx/src/share/tests/combine_ops.cpp

@@ -0,0 +1,74 @@
+#include <catch2/catch.hpp>
+
+#include "share/util/eamxx_combine_ops.hpp"
+#include "share/util/eamxx_universal_constants.hpp"
+#include "share/eamxx_types.hpp"
+
+#include "ekat/ekat_pack.hpp"
+
+#include <iomanip>
+
+namespace {
+
+TEST_CASE ("combine_ops") {
+  using namespace scream;
+  using pack_type = ekat::Pack<Real,SCREAM_PACK_SIZE>;
+
+  constexpr auto Replace   = CombineMode::Replace;
+  constexpr auto Update    = CombineMode::Update;
+  constexpr auto Multiply  = CombineMode::Multiply;
+  constexpr auto Divide    = CombineMode::Divide;
+  constexpr auto Max       = CombineMode::Max;
+  constexpr auto Min       = CombineMode::Min;
+  constexpr auto fv_val    = constants::fill_value<Real>;
+
+  const pack_type two (2.0);
+  const pack_type four (4.0);
+  const pack_type six (6.0);
+  const pack_type ten (10.0);
+  const pack_type fv (constants::fill_value<Real>);
+  const pack_type fv_times_ten (constants::fill_value<Real> * 10);
+
+  pack_type x;
+
+  x = two;
+  combine<Replace>(two,x,1,0);
+  REQUIRE ( (x==two).all() );
+
+  combine<Update>(two,x,2.0,1.0);
+  REQUIRE ( (x==six).all() );
+  fill_aware_combine<Update>(fv,x,fv_val,2.0,1.0);
+  if (not (x==six).all() ) {
+    std::cout << "x: " << x << "\n";
+    std::cout << " x[0]: " << std::setprecision(18) << x[0] << "\n";
+    std::cout << "fv[0]: " << std::setprecision(18) << fv[0] << "\n";
+  }
+  REQUIRE ( (x==six).all() );
+
+  x = two;
+  combine<Multiply>(two,x,1,1);
+  REQUIRE ( (x==four).all() );
+  fill_aware_combine<Multiply>(fv,x,fv_val,1,1);
+  REQUIRE ( (x==four).all() );
+
+  x = four;
+  combine<Divide>(two,x,1,1);
+  REQUIRE ( (x==two).all() );
+  fill_aware_combine<Divide>(fv,x,fv_val,1,1);
+  REQUIRE ( (x==two).all() );
+
+  x = two;
+  combine<Max>(four,x,1,1);
+  REQUIRE ( (x==four).all() );
+  fill_aware_combine<Max>(fv,x,fv_val,1,1);
+  REQUIRE ( (x==four).all() );
+
+  x = four;
+  combine<Min>(two,x,1,1);
+  REQUIRE ( (x==two).all() );
+  x = fv_times_ten;
+  fill_aware_combine<Min>(fv,x,fv_val,1,1);
+  REQUIRE ( (x==fv_times_ten).all() );
+}
+
+} // anonymous namespace

components/eamxx/src/share/tests/field_tests.cpp

@@ -887,6 +887,13 @@ TEST_CASE ("update") {
       Field f4 = two.clone();
       f4.min(f3);
       REQUIRE (views_are_equal(f3, f4));
+
+      // Check that updating with rhs==fill_value ignores the rhs
+      f3.deep_copy(constants::fill_value<Real>);
+      f3.get_header().set_extra_data("mask_value",constants::fill_value<Real>);
+      f2.deep_copy(1.0);
+      f2.max(f3);
+      REQUIRE (views_are_equal(f2,one));
     }
 
     SECTION ("int") {
@@ -904,6 +911,13 @@ TEST_CASE ("update") {
       Field f4 = two.clone();
       f4.min(f3);
       REQUIRE (views_are_equal(f3, f4));
+
+      // Check that updating with rhs==fill_value ignores the rhs
+      f3.deep_copy(constants::fill_value<int>);
+      f3.get_header().set_extra_data("mask_value",constants::fill_value<int>);
+      f2.deep_copy(1);
+      f2.max(f3);
+      REQUIRE (views_are_equal(f2,one));
     }
   }
 
@@ -949,6 +963,19 @@ TEST_CASE ("update") {
       // Same, but we discard current content of f3
       f3.update(f_real,2,0);
       REQUIRE (views_are_equal(f3,f2));
+
+      // Check that updating with rhs==fill_value ignores the rhs
+      Field one = f_real.clone();
+      one.deep_copy(1.0);
+
+      f3.deep_copy(constants::fill_value<Real>);
+      f3.get_header().set_extra_data("mask_value",constants::fill_value<Real>);
+      f2.deep_copy(1.0);
+      f2.update(f3,1,1);
+      if (not views_are_equal(f2,one)) {
+        print_field_hyperslab(f2);
+      }
+      REQUIRE (views_are_equal(f2,one));
     }
 
     SECTION ("int") {
@@ -968,6 +995,16 @@ TEST_CASE ("update") {
       // Same, but we discard current content of f3
       f3.update(f_int,2,0);
       REQUIRE (views_are_equal(f3,f2));
+
+      // Check that updating with rhs==fill_value ignores the rhs
+      Field one = f_int.clone();
+      one.deep_copy(1);
+
+      f3.deep_copy(constants::fill_value<int>);
+      f3.get_header().set_extra_data("mask_value",constants::fill_value<int>);
+      f2.deep_copy(1);
+      f2.update(f3,1,1);
+      REQUIRE (views_are_equal(f2,one));
     }
   }
 }

components/eamxx/src/share/tests/utils_tests.cpp

@@ -7,6 +7,16 @@
 #include "share/util/eamxx_setup_random_test.hpp"
 #include "share/eamxx_config.hpp"
 
+TEST_CASE("fill_value") {
+  using namespace scream::constants;
+
+  // Ensure we have the SAME numerical value for both float and double
+  auto fv_d = fill_value<double>;
+  auto fv_f = fill_value<float>;
+
+  REQUIRE (fv_d==fv_f);
+}
+
 TEST_CASE("contiguous_superset") {
   using namespace scream;
 

components/eamxx/src/share/util/eamxx_combine_ops.hpp

@@ -5,6 +5,7 @@
 
 #include <ekat/ekat_scalar_traits.hpp>
 #include <ekat/util/ekat_math_utils.hpp>
+#include <ekat/util/ekat_where.hpp>
 
 // For KOKKOS_INLINE_FUNCTION
 #include <Kokkos_Core.hpp>
@@ -75,28 +76,39 @@ void combine (const ScalarIn& newVal, ScalarOut& result,
       break;
   }
 }
-/* Special version of combine that takes a mask into account */
+
+// Special version of combine that ignores newVal if newVal==fill_value.
+// Replace is the only combine mode that is allowed to consider fill_val values.
+// This is b/c it's the only way we can use this function inside Field method/utils
+// in order to set all entries of a Field to fill_val. You can also think of fill_val
+// as a special number for which the arithmetic operations are not defined.
+// All CM except for Replace involve an arithmetic op between of two numbers,
+// so combining with fill_val makes no sense. However, it makes sense to set
+// an output variable to fill_val.
 template<CombineMode CM, typename ScalarIn, typename ScalarOut,
          typename CoeffType = typename ekat::ScalarTraits<ScalarIn>::scalar_type>
 KOKKOS_FORCEINLINE_FUNCTION
-void fill_aware_combine (const ScalarIn& newVal, ScalarOut& result, const ScalarOut fill_val,
+void fill_aware_combine (const ScalarIn& newVal, ScalarOut& result,
+              const typename ekat::ScalarTraits<ScalarIn>::scalar_type fill_val,
               const CoeffType alpha, const CoeffType beta)
 {
-  switch (CM) {
-    case CombineMode::Replace:
-      combine<CM>(newVal,result,alpha,beta);
-      break;
-    case CombineMode::Update:
-    case CombineMode::Multiply:
-    case CombineMode::Divide:
-    case CombineMode::Max:
-    case CombineMode::Min:
-      if (newVal != fill_val)
-        combine<CM>(newVal,result,alpha,beta);
-      break;
-
-    default:
-      EKAT_KERNEL_ERROR_MSG("Unsupported combine mode for 'fill_aware_combine' overload");
+  if constexpr (CM==CombineMode::Replace) {
+    return combine<CM>(newVal,result,alpha,beta);
+  }
+
+  // The where object will perform the assignment ONLY where the mask is true
+  auto where = ekat::where(newVal!=fill_val,result);
+  if (where.any()) {
+    // TODO: I thought about doing the switch manually, and do stuff like (e.g., for Update)
+    //  where *= beta;
+    //  where += alpha*newVal
+    // but there is no packed version of where.max(rhs), only a scalar version
+    // (meaning a version where rhs is a scalar, not a pack).
+    // If ekat::where_expression ever implements a packed overload for max/min,
+    // we can get rid of the temporary by doing a manual switch.
+    auto tmp = result;
+    combine<CM>(newVal,tmp,alpha,beta);
+    where = tmp;
   }
 }