diff --git a/components/eamxx/src/diagnostics/aodvis.cpp b/components/eamxx/src/diagnostics/aodvis.cpp index 0c3f7f806d8c..2520d18ac8f1 100644 --- a/components/eamxx/src/diagnostics/aodvis.cpp +++ b/components/eamxx/src/diagnostics/aodvis.cpp @@ -48,9 +48,6 @@ 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::fill_value; - - m_mask_val = m_params.get("mask_value", var_fill_value); std::string mask_name = name() + " mask"; FieldLayout mask_layout({COL}, {m_ncols}); @@ -58,8 +55,7 @@ void AODVis::initialize_impl(const RunType /*run_type*/) { Field diag_mask(mask_fid); diag_mask.allocate_view(); - m_diagnostic_output.get_header().set_extra_data("mask_data", diag_mask); - m_diagnostic_output.get_header().set_extra_data("mask_value", m_mask_val); + m_diagnostic_output.get_header().set_extra_data("mask_field", diag_mask); } void AODVis::compute_diagnostic_impl() { @@ -68,9 +64,11 @@ void AODVis::compute_diagnostic_impl() { using TPF = ekat::TeamPolicyFactory; using RU = ekat::ReductionUtils; + constexpr auto fill_value = constants::fill_value; + const auto aod = m_diagnostic_output.get_view(); const auto mask = m_diagnostic_output.get_header() - .get_extra_data("mask_data") + .get_extra_data("mask_field") .get_view(); const auto tau_vis = get_field_in("aero_tau_sw") .subfield(1, m_vis_bnd) @@ -78,13 +76,12 @@ void AODVis::compute_diagnostic_impl() { const auto sunlit = get_field_in("sunlit").get_view(); const auto num_levs = m_nlevs; - const auto var_fill_value = m_mask_val; const auto policy = TPF::get_default_team_policy(m_ncols, m_nlevs); Kokkos::parallel_for( "Compute " + m_diagnostic_output.name(), policy, KOKKOS_LAMBDA(const MT &team) { const int icol = team.league_rank(); if(sunlit(icol) == 0.0) { - aod(icol) = var_fill_value; + aod(icol) = fill_value; Kokkos::single(Kokkos::PerTeam(team), [&] { mask(icol) = 0; }); } else { auto tau_icol = ekat::subview(tau_vis, icol); diff --git a/components/eamxx/src/diagnostics/aodvis.hpp b/components/eamxx/src/diagnostics/aodvis.hpp index d366d0ed6376..00be05a97f18 100644 --- a/components/eamxx/src/diagnostics/aodvis.hpp +++ b/components/eamxx/src/diagnostics/aodvis.hpp @@ -36,8 +36,6 @@ class AODVis : public AtmosphereDiagnostic { int m_swbands = eamxx_swbands(); int m_vis_bnd = eamxx_vis_swband_idx(); - - Real m_mask_val; }; } // namespace scream diff --git a/components/eamxx/src/diagnostics/atm_backtend.cpp b/components/eamxx/src/diagnostics/atm_backtend.cpp index f31b3c76158d..00f06f83e79b 100644 --- a/components/eamxx/src/diagnostics/atm_backtend.cpp +++ b/components/eamxx/src/diagnostics/atm_backtend.cpp @@ -61,7 +61,6 @@ void AtmBackTendDiag::init_timestep(const util::TimeStamp &start_of_step) { } void AtmBackTendDiag::compute_diagnostic_impl() { - Real var_fill_value = constants::fill_value; std::int64_t dt; const auto &f = get_field_in(m_name); @@ -77,7 +76,7 @@ void AtmBackTendDiag::compute_diagnostic_impl() { } else { // This is the first time we evaluate this diag. We cannot compute a tend // yet, so fill with an invalid value - m_diagnostic_output.deep_copy(var_fill_value); + m_diagnostic_output.deep_copy(constants::fill_value); } } diff --git a/components/eamxx/src/diagnostics/field_at_pressure_level.cpp b/components/eamxx/src/diagnostics/field_at_pressure_level.cpp index d44704ab5341..a86f427d07c5 100644 --- a/components/eamxx/src/diagnostics/field_at_pressure_level.cpp +++ b/components/eamxx/src/diagnostics/field_at_pressure_level.cpp @@ -85,16 +85,14 @@ 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("mask_value",Real(constants::fill_value)); - std::string mask_name = m_diag_name + " mask"; FieldLayout mask_layout( {COL}, {num_cols}); FieldIdentifier mask_fid (mask_name,mask_layout, nondim, gname); Field diag_mask(mask_fid); diag_mask.allocate_view(); - m_diagnostic_output.get_header().set_extra_data("mask_data",diag_mask); - m_diagnostic_output.get_header().set_extra_data("mask_value",m_mask_val); + m_diagnostic_output.get_header().set_extra_data("mask_field",diag_mask); + m_diagnostic_output.get_header().set_may_be_filled(true); using stratts_t = std::map; @@ -126,11 +124,11 @@ void FieldAtPressureLevel::compute_diagnostic_impl() const int nlevs = pl.dim(1); auto p_tgt = m_pressure_level; - auto mval = m_mask_val; + constexpr auto fval = constants::fill_value; if (rank==2) { auto policy = KT::RangePolicy(0,ncols); auto diag = m_diagnostic_output.get_view(); - auto mask = m_diagnostic_output.get_header().get_extra_data("mask_data").get_view(); + auto mask = m_diagnostic_output.get_header().get_extra_data("mask_field").get_view(); auto f_v = f.get_view(); Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int icol) { auto x1 = ekat::subview(p_src_v,icol); @@ -139,7 +137,7 @@ void FieldAtPressureLevel::compute_diagnostic_impl() auto end = beg + nlevs; auto last = beg + (nlevs-1); if (p_tgt<*beg or p_tgt>*last) { - diag(icol) = mval; + diag(icol) = fval; mask(icol) = 0; } else { auto ub = ekat::upper_bound(beg,end,p_tgt); @@ -161,7 +159,7 @@ void FieldAtPressureLevel::compute_diagnostic_impl() const int ndims = f.get_header().get_identifier().get_layout().get_vector_dim(); auto policy = KT::TeamPolicy(ncols,ndims); auto diag = m_diagnostic_output.get_view(); - auto mask = m_diagnostic_output.get_header().get_extra_data("mask_data").get_view(); + auto mask = m_diagnostic_output.get_header().get_extra_data("mask_field").get_view(); auto f_v = f.get_view(); Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const MemberType& team) { int icol = team.league_rank(); @@ -171,7 +169,7 @@ void FieldAtPressureLevel::compute_diagnostic_impl() auto last = beg + (nlevs-1); Kokkos::parallel_for(Kokkos::TeamVectorRange(team,ndims),[&](const int idim) { if (p_tgt<*beg or p_tgt>*last) { - diag(icol,idim) = mval; + diag(icol,idim) = fval; Kokkos::single(Kokkos::PerTeam(team),[&]{ mask(icol) = 0; }); diff --git a/components/eamxx/src/diagnostics/field_at_pressure_level.hpp b/components/eamxx/src/diagnostics/field_at_pressure_level.hpp index 6a6f687b331e..4a0e39e194b9 100644 --- a/components/eamxx/src/diagnostics/field_at_pressure_level.hpp +++ b/components/eamxx/src/diagnostics/field_at_pressure_level.hpp @@ -37,8 +37,6 @@ class FieldAtPressureLevel : public AtmosphereDiagnostic Real m_pressure_level; int m_num_levs; - Real m_mask_val; - }; // class FieldAtPressureLevel } //namespace scream diff --git a/components/eamxx/src/diagnostics/tests/aodvis_test.cpp b/components/eamxx/src/diagnostics/tests/aodvis_test.cpp index bba86648a058..3acb550fe57b 100644 --- a/components/eamxx/src/diagnostics/tests/aodvis_test.cpp +++ b/components/eamxx/src/diagnostics/tests/aodvis_test.cpp @@ -32,8 +32,6 @@ TEST_CASE("aodvis") { using namespace ShortFieldTagsNames; using namespace ekat::units; - Real var_fill_value = constants::fill_value; - Real some_limit = 0.0025; // A world comm @@ -121,7 +119,7 @@ TEST_CASE("aodvis") { for(int icol = 0; icol < grid->get_num_local_dofs(); ++icol) { if(sun_h(icol) < some_limit) { - aod_t(icol) = var_fill_value; + aod_t(icol) = constants::fill_value; } else { for(int ilev = 0; ilev < nlevs; ++ilev) { aod_t(icol) += tau_h(icol, swvis, ilev); diff --git a/components/eamxx/src/diagnostics/tests/atm_backtend_test.cpp b/components/eamxx/src/diagnostics/tests/atm_backtend_test.cpp index 734604aca2d4..914b30ab6139 100644 --- a/components/eamxx/src/diagnostics/tests/atm_backtend_test.cpp +++ b/components/eamxx/src/diagnostics/tests/atm_backtend_test.cpp @@ -65,8 +65,6 @@ TEST_CASE("atm_backtend") { REQUIRE_THROWS(diag_factory.create("AtmBackTendDiag", comm, params)); // No 'tendency_name' - Real var_fill_value = constants::fill_value; - // Set time for qc and randomize its values qc.get_header().get_tracking().update_time_stamp(t0); randomize(qc, engine, pdf); @@ -83,9 +81,9 @@ TEST_CASE("atm_backtend") { diag->compute_diagnostic(); auto diag_f = diag->get_diagnostic(); - // Check result: diag should be filled with var_fill_value + // Check result: diag should be filled with fill_value auto some_field = qc.clone(); - some_field.deep_copy(var_fill_value); + some_field.deep_copy(constants::fill_value); REQUIRE(views_are_equal(diag_f, some_field)); const Real a_day = 24.0 * 60.0 * 60.0; // seconds diff --git a/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp b/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp index 991bcded9a39..1c17ae4d1e4c 100644 --- a/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp +++ b/components/eamxx/src/diagnostics/tests/field_at_pressure_level_tests.cpp @@ -104,7 +104,7 @@ TEST_CASE("field_at_pressure_level_p2") diag_f.sync_to_host(); auto test2_diag_v = diag_f.get_view(); // Check the mask field inside the diag_f - auto mask_f = diag_f.get_header().get_extra_data("mask_data"); + auto mask_f = diag_f.get_header().get_extra_data("mask_field"); mask_f.sync_to_host(); auto test2_mask_v = mask_f.get_view(); // @@ -125,13 +125,12 @@ TEST_CASE("field_at_pressure_level_p2") diag_f.sync_to_host(); auto test2_diag_v = diag_f.get_view(); // Check the mask field inside the diag_f - auto mask_f = diag_f.get_header().get_extra_data("mask_data"); + auto mask_f = diag_f.get_header().get_extra_data("mask_field"); mask_f.sync_to_host(); auto test2_mask_v = mask_f.get_view(); - auto mask_val = diag_f.get_header().get_extra_data("mask_value"); - // + for (int icol=0;icol)); REQUIRE(approx(test2_mask_v(icol),Real(0.0))); } } diff --git a/components/eamxx/src/physics/nudging/eamxx_nudging_process_interface.cpp b/components/eamxx/src/physics/nudging/eamxx_nudging_process_interface.cpp index 93787083f361..fad415d7e672 100644 --- a/components/eamxx/src/physics/nudging/eamxx_nudging_process_interface.cpp +++ b/components/eamxx/src/physics/nudging/eamxx_nudging_process_interface.cpp @@ -344,21 +344,17 @@ void Nudging::run_impl (const double dt) const auto fl = f.get_header().get_identifier().get_layout(); const auto v = f.get_view(); - Real var_fill_value = constants::fill_value; - // 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("mask_value"); - } + constexpr Real fill_value = constants::fill_value; const int ncols = fl.dim(0); const int nlevs = fl.dim(1); - const auto thresh = std::abs(var_fill_value)*0.0001; + const auto thresh = std::abs(fill_value)*0.0001; auto lambda = KOKKOS_LAMBDA(const int icol) { int first_good = nlevs; int last_good = -1; for (int k=0; kthresh) { - // This entry is substantially different from var_fill_value, so it's good + if (std::abs(v(icol,k)-fill_value)>thresh) { + // This entry is substantially different from fill_value, so it's good first_good = ekat::impl::min(first_good,k); last_good = ekat::impl::max(last_good,k); } diff --git a/components/eamxx/src/physics/nudging/tests/nudging_tests_helpers.hpp b/components/eamxx/src/physics/nudging/tests/nudging_tests_helpers.hpp index a2aafb856abc..05c69f77ae7f 100644 --- a/components/eamxx/src/physics/nudging/tests/nudging_tests_helpers.hpp +++ b/components/eamxx/src/physics/nudging/tests/nudging_tests_helpers.hpp @@ -13,7 +13,6 @@ constexpr int nlevs_data = 20; constexpr int nsteps_data = 5; constexpr int dt_data = 100; constexpr int nlevs_filled = 2; -constexpr double fill_val = 1e30; util::TimeStamp get_t0 () { return util::TimeStamp ({2000,1,1},{12,0,0}); @@ -89,27 +88,27 @@ void compute_field (Field f, const auto f_h = f.get_view(); for (int icol=0;icol; } for (int ilev=lev_beg;ilev; } } } else { const auto f_h = f.get_view(); for (int icol=0;icol; } for (int ilev=lev_beg;ilev; } } } @@ -150,7 +149,6 @@ create_om (const std::string& filename_prefix, params.set("filename_prefix",filename_prefix); params.set("floating_point_precision","real"); params.set("field_names",strvec_t{"p_mid","U","V"}); - params.set("fill_value",fill_val); auto& ctrl_pl = params.sublist("output_control"); ctrl_pl.set("frequency_units","nsteps"); diff --git a/components/eamxx/src/share/field/field_header.hpp b/components/eamxx/src/share/field/field_header.hpp index 125b53f72704..3f167f17d5a9 100644 --- a/components/eamxx/src/share/field/field_header.hpp +++ b/components/eamxx/src/share/field/field_header.hpp @@ -87,6 +87,9 @@ class FieldHeader : public FamilyTracking { // - they have the same tracking, alloc_prop and extra data (they were created by alias above) // - they have the same parent field and their subview info (form alloc prop) are the same bool is_aliasing (const FieldHeader& rhs) const; + + bool may_be_filled () const { return has_extra_data("may_be_filled") and get_extra_data("may_be_filled"); } + void set_may_be_filled (const bool value) { set_extra_data("may_be_filled",value); } protected: // Friend this function, so it can set up a subfield header diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index a5365415c872..3db6991e801b 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -559,17 +559,17 @@ update (const Field& x, const ST alpha, const ST beta) " - x layout: " + x_l.to_string() + "\n" " - y layout: " + y_l.to_string() + "\n"); - // Determine if the RHS has fill value entries that require extra treatment - bool use_fill = x.get_header().has_extra_data("mask_value"); + // Determine if the RHS can contain fill_value entries + bool fill_aware = x.get_header().may_be_filled(); if (dt==DataType::IntType) { - if (use_fill) { + if (fill_aware) { return update_impl(x,alpha,beta); } else { return update_impl(x,alpha,beta); } } else if (dt==DataType::FloatType) { - if (use_fill) { + if (fill_aware) { if (rhs_dt==DataType::FloatType) return update_impl(x,alpha,beta); else @@ -581,7 +581,7 @@ update (const Field& x, const ST alpha, const ST beta) return update_impl(x,alpha,beta); } } else if (dt==DataType::DoubleType) { - if (use_fill) { + if (fill_aware) { if (rhs_dt==DataType::DoubleType) return update_impl(x,alpha,beta); else if (rhs_dt==DataType::FloatType) @@ -601,151 +601,142 @@ update (const Field& x, const ST alpha, const ST beta) } } -template +template void Field:: update_impl (const Field& x, const ST alpha, const ST beta) { const auto& layout = x.get_header().get_identifier().get_layout(); const auto& dims = layout.dims(); - XST fill_val = 0; - if constexpr (use_fill) { - EKAT_REQUIRE_MSG (x.get_header().has_extra_data("mask_value"), - "Error! Field::update_impl called with use_fill, but x does not have mask_value extra data.\n" - " - *this name: " + name() + "\n" - " - x name: " + x.name() + "\n"); - fill_val = x.get_header().get_extra_data("mask_value"); - } - // Must handle the case where one of the two views is strided (or both) const auto x_contig = x.get_header().get_alloc_properties().contiguous(); const auto y_contig = get_header().get_alloc_properties().contiguous(); switch (layout.rank()) { case 0: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; case 1: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; case 2: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; case 3: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; case 4: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; case 5: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; case 6: if (x_contig and y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (x_contig) - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else if (y_contig) - details::cvh(get_view(), + details::cvh(get_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); else - details::cvh(get_strided_view(), + details::cvh(get_strided_view(), x.get_strided_view(), - alpha,beta,fill_val,dims); + alpha,beta,dims); break; default: EKAT_ERROR_MSG ("Error! Rank not supported in update_field.\n" diff --git a/components/eamxx/src/share/field/field_impl_details.hpp b/components/eamxx/src/share/field/field_impl_details.hpp index 32ddb7e64447..d0751861e184 100644 --- a/components/eamxx/src/share/field/field_impl_details.hpp +++ b/components/eamxx/src/share/field/field_impl_details.hpp @@ -10,7 +10,7 @@ namespace scream namespace details { -template +template struct CombineViewsHelper { using exec_space = typename LhsView::traits::execution_space; @@ -52,77 +52,54 @@ struct CombineViewsHelper { KOKKOS_INLINE_FUNCTION void operator() (int i) const { - if constexpr (use_fill) - if constexpr (N==0) - fill_aware_combine(rhs(),lhs(),fill_val,alpha,beta); - else - fill_aware_combine(rhs(i),lhs(i),fill_val,alpha,beta); + if constexpr (N==0) + combine(rhs(),lhs(),alpha,beta); else - if constexpr (N==0) - combine(rhs(),lhs(),alpha,beta); - else - combine(rhs(i),lhs(i),alpha,beta); + combine(rhs(i),lhs(i),alpha,beta); } KOKKOS_INLINE_FUNCTION void operator() (int i, int j) const { - if constexpr (use_fill) - fill_aware_combine(rhs(i,j),lhs(i,j),fill_val,alpha,beta); - else - combine(rhs(i,j),lhs(i,j),alpha,beta); + combine(rhs(i,j),lhs(i,j),alpha,beta); } KOKKOS_INLINE_FUNCTION void operator() (int i, int j, int k) const { - if constexpr (use_fill) - fill_aware_combine(rhs(i,j,k),lhs(i,j,k),fill_val,alpha,beta); - else - combine(rhs(i,j,k),lhs(i,j,k),alpha,beta); + combine(rhs(i,j,k),lhs(i,j,k),alpha,beta); } KOKKOS_INLINE_FUNCTION void operator() (int i, int j, int k, int l) const { - if constexpr (use_fill) - fill_aware_combine(rhs(i,j,k,l),lhs(i,j,k,l),fill_val,alpha,beta); - else - combine(rhs(i,j,k,l),lhs(i,j,k,l),alpha,beta); + combine(rhs(i,j,k,l),lhs(i,j,k,l),alpha,beta); } KOKKOS_INLINE_FUNCTION void operator() (int i, int j, int k, int l, int m) const { - if constexpr (use_fill) - fill_aware_combine(rhs(i,j,k,l,m),lhs(i,j,k,l,m),fill_val,alpha,beta); - else - combine(rhs(i,j,k,l,m),lhs(i,j,k,l,m),alpha,beta); + combine(rhs(i,j,k,l,m),lhs(i,j,k,l,m),alpha,beta); } KOKKOS_INLINE_FUNCTION void operator() (int i, int j, int k, int l, int m, int n) const { - if constexpr (use_fill) - fill_aware_combine(rhs(i,j,k,l,m,n),lhs(i,j,k,l,m,n),fill_val,alpha,beta); - else - combine(rhs(i,j,k,l,m,n),lhs(i,j,k,l,m,n),alpha,beta); + combine(rhs(i,j,k,l,m,n),lhs(i,j,k,l,m,n),alpha,beta); } ST alpha; ST beta; LhsView lhs; RhsView rhs; - typename LhsView::traits::value_type fill_val; }; -template +template void cvh (LhsView lhs, RhsView rhs, - ST alpha, ST beta, typename RhsView::traits::value_type fill_val, + ST alpha, ST beta, const std::vector& dims) { - CombineViewsHelper helper; + CombineViewsHelper helper; helper.lhs = lhs; helper.rhs = rhs; helper.alpha = alpha; helper.beta = beta; - helper.fill_val = fill_val; helper.run(dims); } diff --git a/components/eamxx/src/share/grid/remap/coarsening_remapper.cpp b/components/eamxx/src/share/grid/remap/coarsening_remapper.cpp index c138d872b07b..b4ceadbceb9b 100644 --- a/components/eamxx/src/share/grid/remap/coarsening_remapper.cpp +++ b/components/eamxx/src/share/grid/remap/coarsening_remapper.cpp @@ -96,18 +96,13 @@ registration_ends_impl () for (int i=0; i("mask_data"); - const auto& src_mask_val = src.get_header().get_extra_data("mask_value"); + const auto& src_mask = src.get_header().get_extra_data("mask_field"); // Make sure fields representing masks are not themselves meant to be masked. - EKAT_REQUIRE_MSG(not src_mask.get_header().has_extra_data("mask_data"), + EKAT_REQUIRE_MSG(not src_mask.get_header().has_extra_data("mask_field"), "Error! A mask field cannot be itself masked.\n" " - field name: " + src.name() + "\n" " - mask field name: " + src_mask.name() + "\n"); @@ -127,20 +122,6 @@ registration_ends_impl () " - field layout: " + f_lt.to_string() + "\n" " - mask layout: " + m_lt.to_string() + "\n"); - // If not there, set mask value in the tgt field too - if (tgt.get_header().has_extra_data("mask_value")) { - const auto& tgt_mask_val = tgt.get_header().get_extra_data("mask_value"); - - EKAT_REQUIRE_MSG (tgt_mask_val==src_mask_val, - "Error! Target field stores a mask data different from the src field.\n" - " - src field name: " + src.name() + "\n" - " - tgt field name: " + tgt.name() + "\n" - " - src mask value: " << src_mask_val << "\n" - " - tgt mask value: " << tgt_mask_val << "\n"); - } else { - tgt.get_header().set_extra_data("mask_value",src_mask_val); - } - // If it's the first time we find this mask, store it, so we can register later const auto& src_mask_fid = src_mask.get_header().get_identifier(); int mask_idx = get_mask_idx(src_mask_fid); @@ -153,7 +134,7 @@ registration_ends_impl () masks.push_back(std::make_pair(src_mask,tgt_mask)); mask_idx = masks.size()-1; } - tgt.get_header().set_extra_data("mask_data",masks[mask_idx].second); + tgt.get_header().set_extra_data("mask_field",masks[mask_idx].second); } // Add all masks to the fields to remap @@ -196,10 +177,10 @@ void CoarseningRemapper::remap_fwd_impl () const auto& f_src = m_src_fields[i]; const auto& f_ov = m_ov_fields[i]; - const bool masked = m_track_mask and f_src.get_header().has_extra_data("mask_data"); + const bool masked = m_track_mask and f_src.get_header().has_extra_data("mask_field"); if (masked) { // Pass the mask to the local_mat_vec routine - const auto& mask = f_src.get_header().get_extra_data("mask_data"); + const auto& mask = f_src.get_header().get_extra_data("mask_field"); // If possible, dispatch kernel with SCREAM_PACK_SIZE if (can_pack_field(f_src) and can_pack_field(f_ov) and can_pack_field(mask)) { @@ -235,9 +216,9 @@ void CoarseningRemapper::remap_fwd_impl () if (m_track_mask) { for (int i=0; i("mask_data"); + const auto& mask = f_tgt.get_header().get_extra_data("mask_field"); if (can_pack_field(f_tgt) and can_pack_field(mask)) { rescale_masked_fields(f_tgt,mask); } else { @@ -258,15 +239,11 @@ rescale_masked_fields (const Field& x, const Field& mask) const using Pack = ekat::Pack; using PackInfo = ekat::PackInfo; + constexpr auto fill_val = constants::fill_value; + const auto& layout = x.get_header().get_identifier().get_layout(); const int rank = layout.rank(); const int ncols = m_tgt_grid->get_num_local_dofs(); - Real mask_val = std::numeric_limits::max()/10.0; - if (x.get_header().has_extra_data("mask_value")) { - mask_val = x.get_header().get_extra_data("mask_value"); - } else { - EKAT_ERROR_MSG ("ERROR! Field " + x.name() + " is masked, but stores no mask_value extra data.\n"); - } const Real mask_threshold = std::numeric_limits::epsilon(); // TODO: Should we not hardcode the threshold for simply masking out the column. switch (rank) { @@ -282,7 +259,7 @@ rescale_masked_fields (const Field& x, const Field& mask) const if (m_view(icol)>mask_threshold) { x_view(icol) /= m_view(icol); } else { - x_view(icol) = mask_val; + x_view(icol) = fill_val; } }); break; @@ -322,7 +299,7 @@ rescale_masked_fields (const Field& x, const Field& mask) const if (masked.any()) { x_sub(j).set(masked,x_sub(j)/m_sub(j)); } - x_sub(j).set(!masked,mask_val); + x_sub(j).set(!masked,fill_val); }); } }); @@ -370,7 +347,7 @@ rescale_masked_fields (const Field& x, const Field& mask) const if (masked.any()) { x_sub(k).set(masked,x_sub(k)/m_sub(k)); } - x_sub(k).set(!masked,mask_val); + x_sub(k).set(!masked,fill_val); }); } }); @@ -421,7 +398,7 @@ rescale_masked_fields (const Field& x, const Field& mask) const if (masked.any()) { x_sub(l).set(masked,x_sub(l)/m_sub(l)); } - x_sub(l).set(!masked,mask_val); + x_sub(l).set(!masked,fill_val); }); } }); diff --git a/components/eamxx/src/share/grid/remap/vertical_remapper.cpp b/components/eamxx/src/share/grid/remap/vertical_remapper.cpp index 8d5305f2a6c8..4726c26ec337 100644 --- a/components/eamxx/src/share/grid/remap/vertical_remapper.cpp +++ b/components/eamxx/src/share/grid/remap/vertical_remapper.cpp @@ -84,15 +84,6 @@ set_extrapolation_type (const ExtrapType etype, const TopBot where) } } -void VerticalRemapper:: -set_mask_value (const Real mask_val) -{ - EKAT_REQUIRE_MSG (not Kokkos::isnan(mask_val), - "[VerticalRemapper::set_mask_value] Error! Input mask value must be a valid number.\n"); - - m_mask_val = mask_val; -} - void VerticalRemapper:: set_source_pressure (const Field& p, const ProfileType ptype) { @@ -169,8 +160,7 @@ registration_ends_impl () const auto& src = m_src_fields[i]; auto& tgt = m_tgt_fields[i]; - // Clone src layout, since we may strip dims later for mask creation - auto src_layout = src.get_header().get_identifier().get_layout().clone(); + const auto& src_layout = src.get_header().get_identifier().get_layout().clone(); if (src_layout.has_tag(LEV) or src_layout.has_tag(ILEV)) { // Determine if this field can be handled with packs, and whether it's at midpoints @@ -194,69 +184,47 @@ registration_ends_impl () // NOTE: for now we assume that masking is determined only by the COL,LEV location in space // and that fields with multiple components will have the same masking for each component // at a specific COL,LEV - src_layout.strip_dims({CMP}); + + auto src_layout_no_cmp = src_layout.clone(); + src_layout_no_cmp.strip_dims({CMP}); + auto tgt_layout = create_tgt_layout(src_layout_no_cmp); // I this mask has already been created, retrieve it, otherwise create it - const auto mask_name = m_tgt_grid->name() + "_" + ekat::join(src_layout.names(),"_") + "_mask"; - Field tgt_mask; - if (m_field2type.count(mask_name)==0) { + // CAVEAT: the tgt layout ALWAYS has LEV as vertical dim tag. But we NEED different masks for + // src fields defined at LEV and ILEV. So use src_layout_no_cmp to craft the mask name + const auto mask_name = m_tgt_grid->name() + "_" + ekat::join(src_layout_no_cmp.names(),"_") + "_mask"; + auto& mask = m_masks[mask_name]; + if (not mask.is_allocated()) { auto nondim = ekat::units::Units::nondimensional(); // Create this src/tgt mask fields, and assign them to these src/tgt fields extra data - FieldIdentifier src_mask_fid (mask_name, src_layout, nondim, m_src_grid->name() ); - FieldIdentifier tgt_mask_fid = create_tgt_fid(src_mask_fid); - - Field src_mask (src_mask_fid); - src_mask.allocate_view(); - - tgt_mask = Field (tgt_mask_fid); - tgt_mask.allocate_view(); - - // Initialize the src mask values to 1.0 - src_mask.deep_copy(1.0); - - m_src_masks.push_back(src_mask); - m_tgt_masks.push_back(tgt_mask); - - auto& mt = m_field2type[src_mask_fid.name()]; - mt.packed = false; - mt.midpoints = src_layout.has_tag(LEV); - } else { - for (size_t i=0; iname() ); + mask = Field (mask_fid); + mask.allocate_view(); } - EKAT_REQUIRE_MSG(not tgt.get_header().has_extra_data("mask_data"), + EKAT_REQUIRE_MSG(not tgt.get_header().has_extra_data("mask_field"), "[VerticalRemapper::registration_ends_impl] Error! Target field already has mask data assigned.\n" " - tgt field name: " + tgt.name() + "\n"); - EKAT_REQUIRE_MSG(not tgt.get_header().has_extra_data("mask_value"), - "[VerticalRemapper::registration_ends_impl] Error! Target field already has mask value assigned.\n" - " - tgt field name: " + tgt.name() + "\n"); - tgt.get_header().set_extra_data("mask_data",tgt_mask); - tgt.get_header().set_extra_data("mask_value",m_mask_val); + tgt.get_header().set_extra_data("mask_field",mask); + + // Since we do mask (at top and/or bot), the tgt field MAY be contain fill_value entries + tgt.get_header().set_may_be_filled(true); } } else { - // If a field does not have LEV or ILEV it may still have mask tracking assigned from somewhere else. + // If a field does not have LEV or ILEV it may still have fill_value tracking assigned from somewhere else. // For instance, this could be a 2d field computed by FieldAtPressureLevel diagnostic. - // In those cases we want to copy that mask tracking to the target field. - if (src.get_header().has_extra_data("mask_data")) { - EKAT_REQUIRE_MSG(not tgt.get_header().has_extra_data("mask_data"), + // In those cases we want to copy that fill_value tracking to the target field. + if (src.get_header().has_extra_data("mask_field")) { + EKAT_REQUIRE_MSG(not tgt.get_header().has_extra_data("mask_field"), "[VerticalRemapper::registration_ends_impl] Error! Target field already has mask data assigned.\n" " - tgt field name: " + tgt.name() + "\n"); - auto src_mask = src.get_header().get_extra_data("mask_data"); - tgt.get_header().set_extra_data("mask_data",src_mask); + auto src_mask = src.get_header().get_extra_data("mask_field"); + tgt.get_header().set_extra_data("mask_field",src_mask); } - if (src.get_header().has_extra_data("mask_value")) { - EKAT_REQUIRE_MSG(not tgt.get_header().has_extra_data("mask_value"), - "[VerticalRemapper::registration_ends_impl] Error! Target field already has mask value assigned.\n" - " - tgt field name: " + tgt.name() + "\n"); - auto src_mask_val = src.get_header().get_extra_data("mask_value"); - tgt.get_header().set_extra_data("mask_value",src_mask_val); + if (src.get_header().may_be_filled()) { + tgt.get_header().set_may_be_filled(true); } } } @@ -397,6 +365,8 @@ create_layout (const FieldLayout& from_layout, void VerticalRemapper::remap_fwd_impl () { + using namespace ShortFieldTagsNames; + // 1. Setup any interp object that was created (if nullptr, no fields need it) if (m_lin_interp_mid_packed) { setup_lin_interp(*m_lin_interp_mid_packed,m_src_pmid,m_tgt_pmid); @@ -411,9 +381,12 @@ void VerticalRemapper::remap_fwd_impl () setup_lin_interp(*m_lin_interp_int_scalar,m_src_pint,m_tgt_pint); } - using namespace ShortFieldTagsNames; + // 2. Init all masks fields (if any) to 1 (signaling no masked entries) + for (auto& [name, mask] : m_masks) { + mask.deep_copy(1); + } - // 2. Interpolate the fields + // 3. Interpolate the fields for (int i=0; i("mask_data"); - auto f_src_mask = f_src.get_header().get_extra_data("mask_data"); + if (f_tgt.get_header().has_extra_data("mask_field")) { + auto f_tgt_mask = f_tgt.get_header().get_extra_data("mask_field"); + auto f_src_mask = f_src.get_header().get_extra_data("mask_field"); f_tgt_mask.deep_copy(f_src_mask); } } } - // 3. Interpolate the mask fields - for (unsigned i=0; i @@ -621,14 +571,15 @@ void VerticalRemapper:: extrapolate (const Field& f_src, const Field& f_tgt, const Field& p_src, - const Field& p_tgt, - const Real mask_val) const + const Field& p_tgt) const { using TPF = ekat::TeamPolicyFactory; using view2d = typename KokkosTypes::view; using view1d = typename KokkosTypes::view; + constexpr auto fill_val = constants::fill_value; + auto src1d = p_src.rank()==1; auto tgt1d = p_tgt.rank()==1; @@ -654,6 +605,12 @@ extrapolate (const Field& f_src, auto etop = m_etype_top; auto ebot = m_etype_bot; auto mid = nlevs_tgt / 2; + auto do_mask = etop==Mask or ebot==Mask; + decltype(f_tgt.get_view()) mask_v; + if (do_mask) { + mask_v = f_tgt.get_header().get_extra_data("mask_field").get_view(); + } + switch(f_src.rank()) { case 2: { @@ -686,7 +643,8 @@ extrapolate (const Field& f_src, if (ebot==P0) { y_tgt[ilev] = y_src[nlevs_src-1]; } else { - y_tgt[ilev] = mask_val; + y_tgt[ilev] = fill_val; + mask_v(icol,ilev) = 0; } } } else { @@ -695,7 +653,8 @@ extrapolate (const Field& f_src, if (etop==P0) { y_tgt[ilev] = y_src[0]; } else { - y_tgt[ilev] = mask_val; + y_tgt[ilev] = fill_val; + mask_v(icol,ilev) = 0; } } } @@ -737,7 +696,8 @@ extrapolate (const Field& f_src, if (ebot==P0) { y_tgt[ilev] = y_src[nlevs_src-1]; } else { - y_tgt[ilev] = mask_val; + y_tgt[ilev] = fill_val; + mask_v(icol,ilev) = 0; } } } else { @@ -746,7 +706,8 @@ extrapolate (const Field& f_src, if (etop==P0) { y_tgt[ilev] = y_src[0]; } else { - y_tgt[ilev] = mask_val; + y_tgt[ilev] = fill_val; + mask_v(icol,ilev) = 0; } } } diff --git a/components/eamxx/src/share/grid/remap/vertical_remapper.hpp b/components/eamxx/src/share/grid/remap/vertical_remapper.hpp index 9f719dfd60e3..264f719d7558 100644 --- a/components/eamxx/src/share/grid/remap/vertical_remapper.hpp +++ b/components/eamxx/src/share/grid/remap/vertical_remapper.hpp @@ -47,7 +47,6 @@ class VerticalRemapper : public AbstractRemapper ~VerticalRemapper () = default; void set_extrapolation_type (const ExtrapType etype, const TopBot where = TopAndBot); - void set_mask_value (const Real mask_val); void set_source_pressure (const Field& p, const ProfileType ptype); void set_target_pressure (const Field& p, const ProfileType ptype); @@ -95,8 +94,7 @@ class VerticalRemapper : public AbstractRemapper const Field& p_src, const Field& p_tgt) const; void extrapolate (const Field& f_src, const Field& f_tgt, - const Field& p_src, const Field& p_tgt, - const Real mask_val) const; + const Field& p_src, const Field& p_tgt) const; template void setup_lin_interp (const ekat::LinInterp& lin_interp, @@ -114,9 +112,8 @@ class VerticalRemapper : public AbstractRemapper ekat::Comm m_comm; - // Source and target fields - std::vector m_src_masks; - std::vector m_tgt_masks; + // Tgt grid masks (in case extrap type at top or bot is Mask) + std::map m_masks; // Vertical profile fields, both for source and target Field m_src_pmid; @@ -135,7 +132,6 @@ class VerticalRemapper : public AbstractRemapper // Extrapolation settings at top/bottom. Default to P0 extrapolation ExtrapType m_etype_top = P0; ExtrapType m_etype_bot = P0; - Real m_mask_val = std::numeric_limits::quiet_NaN(); // We need to remap mid/int fields separately, and we want to use packs if possible, // so we need to divide input fields into 4 separate categories diff --git a/components/eamxx/src/share/io/scorpio_output.cpp b/components/eamxx/src/share/io/scorpio_output.cpp index 631270acc195..6762643c701b 100644 --- a/components/eamxx/src/share/io/scorpio_output.cpp +++ b/components/eamxx/src/share/io/scorpio_output.cpp @@ -13,9 +13,11 @@ #include namespace { - // Helper lambda, to copy io string attributes. This will be used if any - // remapper is created, to ensure atts set by atm_procs are not lost - void transfer_io_str_atts (const scream::Field& src, scream::Field& tgt) { + // Helper lambda, to copy extra data (io string attributes plus filled settings). + // This will be used if any remapper is created, to ensure atts set by atm_procs are not lost + void transfer_extra_data (const scream::Field& src, scream::Field& tgt) { + + // Transfer io string attributes const std::string io_string_atts_key ="io: string attributes"; using stratts_t = std::map; const auto& src_atts = src.get_header().get_extra_data(io_string_atts_key); @@ -23,6 +25,12 @@ namespace { for (const auto& [name,val] : src_atts) { dst_atts[name] = val; } + + // Transfer whether or not this field MAY contain fill_value (to trigger usage of the + // proper implementation of Field update methods + if (src.get_header().may_be_filled()) { + tgt.get_header().set_may_be_filled(true); + } }; } @@ -191,10 +199,6 @@ AtmosphereOutput (const ekat::Comm& comm, const ekat::ParameterList& params, } } - if (params.isParameter("fill_value")) { - m_fill_value = static_cast(params.get("fill_value")); - } - // Setup remappers - if needed auto grid_after_vr = fm_grid; if (use_vertical_remap_from_file) { @@ -204,7 +208,6 @@ AtmosphereOutput (const ekat::Comm& comm, const ekat::ParameterList& params, auto p_int = fm_model->get_field("p_int"); auto vert_remapper = std::make_shared(fm_model->get_grid(),vert_remap_file); vert_remapper->set_source_pressure (p_mid,p_int); - vert_remapper->set_mask_value(m_fill_value); vert_remapper->set_extrapolation_type(VerticalRemapper::Mask); // both Top AND Bot m_vert_remapper = vert_remapper; @@ -214,7 +217,7 @@ AtmosphereOutput (const ekat::Comm& comm, const ekat::ParameterList& params, for (const auto& fname : m_fields_names) { auto src = fm_model->get_field(fname,fm_grid->name()); auto tgt = m_vert_remapper->register_field_from_src(src); - transfer_io_str_atts (src,tgt); + transfer_extra_data (src,tgt); fm_after_vr->add_field(tgt); } m_vert_remapper->registration_ends(); @@ -243,7 +246,7 @@ AtmosphereOutput (const ekat::Comm& comm, const ekat::ParameterList& params, for (const auto& fname : m_fields_names) { auto src = fm_after_vr->get_field(fname,grid_after_vr->name()); auto tgt = m_horiz_remapper->register_field_from_src(src); - transfer_io_str_atts (src,tgt); + transfer_extra_data (src,tgt); fm_after_hr->add_field(tgt); } m_horiz_remapper->registration_ends(); @@ -296,17 +299,14 @@ void AtmosphereOutput::init() // Check if the field for scorpio can alias the field after hremap. // It can do so only for Instant output, and if the field is NOT a subfield ant NOT padded - // Also, if we track avg cnt, we MUST add the mask_value extra data, to trigger fill-value logic + // Also, if we track avg cnt, we MUST add the fill_value extra data, to trigger fill-value logic // when calling Field's update methods if (m_avg_type!=OutputAvgType::Instant or fh.get_alloc_properties().get_padding()>0 or fh.get_parent()!=nullptr) { Field copy(fid); copy.allocate_view(); - transfer_io_str_atts (f,copy); - if (m_track_avg_cnt) { - copy.get_header().set_extra_data("mask_value",Real(m_fill_value)); - } + transfer_extra_data (f,copy); fm_scorpio->add_field(copy); } else { fm_scorpio->add_field(f); @@ -440,11 +440,30 @@ run (const std::string& filename, continue; } + //////////////////////// TODO //////////////////////// + // 1. Make the diags/procs COMPUTE the mask for the field, and + // set the mask as extra data. We DON'T want to compute it here + // 2. Create count with same layout as the mask + // 3. In avg=tally/count, count MAY have smaller layout, hence + // be incompatible. E.g., mask came from FieldAtPressureLevel, + // and was (ncol), but the field is (ncol,2). In this case, we + // ASSUME same mask for all slices, so do scaling on all slices: + // avg.component(i).scale_inv(count) + // 4. In FieldAtPressureLevel, make ALL instances at same Plev share + // the same mask field. How? Two ideas: + // - store a static map string->Field in class, and use a string + // key that encodes grid and press level. Only compute mask if + // timestamp of mask field is "old" + // - make another diag that computes the mask, make F@plev depend + // on that diag. + // 5. FieldAtPressureLevel (and vremap) should only fill the mask field + // if we later need it. E.g, if no AvgCount AND no hremap, we don't need it. + ////////////////////////////////////////////////////// auto field = fm_after_hr->get_field(fname); auto mask = count.get_header().get_extra_data("mask"); - // Find where the field is NOT equal to m_fill_value - compute_mask(field,m_fill_value,mask); + // Find where the field is NOT equal to fill_value + compute_mask(field,constants::fill_value,mask); // mask=1 for "good" entries, and mask=0 otherwise. count.update(mask,1,1); @@ -512,7 +531,7 @@ run (const std::string& filename, f_out.scale_inv(avg_count); const auto& mask = avg_count.get_header().get_extra_data("mask"); - f_out.deep_copy(m_fill_value,mask); + f_out.deep_copy(constants::fill_value,mask); } else { // Divide by steps count only when the summation is complete f_out.scale(Real(1.0) / nsteps_since_last_output); @@ -702,11 +721,9 @@ register_variables(const std::string& filename, // FillValue is a protected metadata, do not add it if it already existed if (fp_precision=="double" or (fp_precision=="real" and std::is_same::value)) { - double fill_value = m_fill_value; - scorpio::set_attribute(filename, alias_name, "_FillValue",fill_value); + scorpio::set_attribute(filename, alias_name, "_FillValue",constants::fill_value); } else { - float fill_value = m_fill_value; - scorpio::set_attribute(filename, alias_name, "_FillValue",fill_value); + scorpio::set_attribute(filename, alias_name, "_FillValue",constants::fill_value); } // If this is has subfields, add list of its children @@ -897,11 +914,11 @@ compute_diagnostics(const bool allow_invalid_fields) if (not computed) { // The diag was either not computable or it may have failed to compute // (e.g., t=0 output with a flux-like diag). - // If we're allowing invalid fields, then we should simply set diag=m_fill_value + // If we're allowing invalid fields, then we should simply set diag=fill_value EKAT_REQUIRE_MSG (allow_invalid_fields, "Error! Failed to compute diagnostic.\n" " - diag name: " + diag->get_diagnostic().name() + "\n"); - d.deep_copy(m_fill_value); + d.deep_copy(constants::fill_value); } } } @@ -950,7 +967,6 @@ init_diagnostics () std::string diag_avg_cnt_name = ""; auto& params = diag->get_params(); if (diag->name()=="FieldAtPressureLevel") { - params.set("mask_value",m_fill_value); diag_avg_cnt_name = "_" + params.get("pressure_value") + params.get("pressure_units"); @@ -963,12 +979,10 @@ init_diagnostics () m_track_avg_cnt |= m_avg_type!=OutputAvgType::Instant; } } else if (diag->name()=="AerosolOpticalDepth550nm") { - params.set("mask_value", m_fill_value); m_track_avg_cnt = m_track_avg_cnt || m_avg_type!=OutputAvgType::Instant; diag_avg_cnt_name = "_" + diag->name(); } else if (diag_field.get_header().has_extra_data("mask_data")) { - params.set("mask_value", m_fill_value); m_track_avg_cnt = m_track_avg_cnt || m_avg_type!=OutputAvgType::Instant; diag_avg_cnt_name = "_" + diag_field.name(); } diff --git a/components/eamxx/src/share/io/scorpio_output.hpp b/components/eamxx/src/share/io/scorpio_output.hpp index f68ff5213b7b..a26875a0042c 100644 --- a/components/eamxx/src/share/io/scorpio_output.hpp +++ b/components/eamxx/src/share/io/scorpio_output.hpp @@ -219,13 +219,6 @@ class AtmosphereOutput DefaultMetadata m_default_metadata; - // Use float, so that if output fp_precision=float, this is a representable value. - // Otherwise, you would get an error from Netcdf, like - // 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::fill_value; - bool m_add_time_dim; bool m_track_avg_cnt = false; std::string m_decomp_dimname = ""; diff --git a/components/eamxx/src/share/io/tests/io_filled.cpp b/components/eamxx/src/share/io/tests/io_filled.cpp index 435d5df8c8a5..2a51faaa4872 100644 --- a/components/eamxx/src/share/io/tests/io_filled.cpp +++ b/components/eamxx/src/share/io/tests/io_filled.cpp @@ -26,7 +26,7 @@ namespace scream { constexpr int num_output_steps = 5; -constexpr Real FillValue = constants::fill_value; +constexpr Real fill_value = constants::fill_value; constexpr Real fill_threshold = 0.5; void set (const Field& f, const double v) { @@ -100,7 +100,7 @@ get_fm (const std::shared_ptr& grid, f.allocate_view(); f.deep_copy(0.0); // For the "filled" field we start with a filled value. f.get_header().get_tracking().update_time_stamp(t0); - f.get_header().set_extra_data("mask_value",FillValue); + f.get_header().set_may_be_filled(true); fm->add_field(f); } @@ -131,7 +131,6 @@ void write (const std::string& avg_type, const std::string& freq_units, om_pl.set("filename_prefix",std::string("io_filled")); om_pl.set("field_names",fnames); om_pl.set("averaging_type", avg_type); - om_pl.set("fill_value",FillValue); om_pl.set("fill_threshold",fill_threshold); om_pl.set("track_avg_cnt",true); auto& ctrl_pl = om_pl.sublist("output_control"); @@ -152,10 +151,10 @@ void write (const std::string& avg_type, const std::string& freq_units, // Update time t += dt; - // Set fields to n+1 or the FillValue, depending on step: + // Set fields to n+1 or the fill_value, depending on step: // - n+1 if n+1 is odd - // - FillValue if n+1 is even - Real setval = ((n+1) % 2 == 0) ? 1.0*(n+1) : FillValue; + // - fill_value if n+1 is even + Real setval = ((n+1) % 2 == 0) ? 1.0*(n+1) : fill_value; for (const auto& n : fnames) { auto f = fm->get_field(n); set(f,setval); @@ -206,7 +205,7 @@ void read (const std::string& avg_type, const std::string& freq_units, reader_pl.set("field_names",fnames); AtmosphereInput reader(reader_pl,fm); - // We set the value n to each input field for each odd valued timestep and FillValue for each even valued timestep + // We set the value n to each input field for each odd valued timestep and fill_value for each even valued timestep // Hence, at output step N = snap*freq, we should get // avg=INSTANT: output = N if (N%2=0), else Fillvalue // avg=MAX: output = N if (N%2=0), else N-1 @@ -232,7 +231,7 @@ void read (const std::string& avg_type, const std::string& freq_units, set(f0,test_val); REQUIRE (views_are_equal(f,f0)); } else if (avg_type=="INSTANT") { - Real test_val = (n*freq%2==0) ? n*freq : FillValue; + Real test_val = (n*freq%2==0) ? n*freq : fill_value; set(f0,test_val); REQUIRE (views_are_equal(f,f0)); } else { // Is avg_type = AVERAGE @@ -243,7 +242,7 @@ void read (const std::string& avg_type, const std::string& freq_units, Real test_val; Real M = freq/2 + (n%2==0 ? 0.0 : 1.0); Real a = n*freq + (n%2==0 ? 0.0 : -1.0); - test_val = (M/freq > fill_threshold) ? a + (M+1.0) : FillValue; + test_val = (M/freq > fill_threshold) ? a + (M+1.0) : fill_value; set(f0,test_val); REQUIRE (views_are_equal(f,f0)); } diff --git a/components/eamxx/src/share/io/tests/io_remap_test.cpp b/components/eamxx/src/share/io/tests/io_remap_test.cpp index b043df9b2f0d..322beada4790 100644 --- a/components/eamxx/src/share/io/tests/io_remap_test.cpp +++ b/components/eamxx/src/share/io/tests/io_remap_test.cpp @@ -16,6 +16,7 @@ namespace { using namespace scream; constexpr int packsize = SCREAM_SMALL_PACK_SIZE; +constexpr Real fill_val = constants::fill_value; using Pack = ekat::Pack; using stratts_t = std::map; @@ -284,8 +285,6 @@ TEST_CASE("io_remap_test","io_remap_test") // --- Vertical Remapping --- { // 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("Fill Value") : constants::fill_value; 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"); @@ -332,7 +331,7 @@ TEST_CASE("io_remap_test","io_remap_test") for (int ii=0; iipi_v(ii,nlevs_src) || p_refpm_v(ii,nlevs_src-1) || p_jjpi_v(ii,nlevs_src) || p_jj("Fill Value") : constants::fill_value; 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"); @@ -431,7 +428,7 @@ TEST_CASE("io_remap_test","io_remap_test") if (found) { Ys_exp /= Ys_wgt; } else { - Ys_exp = mask_val; + Ys_exp = fill_val; } REQUIRE(approx(Ys_v_horiz(ii), Ys_exp)); } @@ -440,8 +437,6 @@ 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("Fill Value") : constants::fill_value; 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"); @@ -502,13 +497,13 @@ TEST_CASE("io_remap_test","io_remap_test") test_mid = (mid_mask_1*calculate_output(p_jj,col1,0)*wgt + mid_mask_2*calculate_output(p_jj,col2,0)*(1-wgt))/(mid_mask_1*wgt + mid_mask_2*(1-wgt)); } else { // This point is completely masked out, assign masked value - test_mid = mask_val; + test_mid = fill_val; } if (int_mask_1 + int_mask_2 > 0.0) { test_int = (int_mask_1*calculate_output(p_jj,col1,0)*wgt + int_mask_2*calculate_output(p_jj,col2,0)*(1-wgt))/(int_mask_1*wgt + int_mask_2*(1-wgt)); } else { // This point is completely masked out, assign masked value - test_int = mask_val; + test_int = fill_val; } REQUIRE(approx(Ym_v_vh(ii,jj), test_mid)); REQUIRE(approx(Yi_v_vh(ii,jj), test_int)); @@ -517,13 +512,13 @@ TEST_CASE("io_remap_test","io_remap_test") test_mid = (mid_mask_1*calculate_output(p_jj,col1,cc+1)*wgt + mid_mask_2*calculate_output(p_jj,col2,cc+1)*(1-wgt))/(mid_mask_1*wgt + mid_mask_2*(1-wgt)); } else { // This point is completely masked out, assign masked value - test_mid = mask_val; + test_mid = fill_val; } if (int_mask_1 + int_mask_2 > 0.0) { test_int = (int_mask_1*calculate_output(p_jj,col1,cc+1)*wgt + int_mask_2*calculate_output(p_jj,col2,cc+1)*(1-wgt))/(int_mask_1*wgt + int_mask_2*(1-wgt)); } else { // This point is completely masked out, assign masked value - test_int = mask_val; + test_int = fill_val; } REQUIRE(approx(Vm_v_vh(ii,cc,jj), test_mid)); REQUIRE(approx(Vi_v_vh(ii,cc,jj), test_int)); @@ -546,7 +541,7 @@ TEST_CASE("io_remap_test","io_remap_test") if (found) { Ys_exp /= Ys_wgt; } else { - Ys_exp = mask_val; + Ys_exp = fill_val; } REQUIRE(approx(Ys_v_vh(ii), Ys_exp)); } @@ -680,7 +675,6 @@ ekat::ParameterList set_output_params(const std::string& name, const std::string params.set("filename_prefix",name); params.set("averaging_type","instant"); - params.set("max_snapshots_per_file",1); params.set("floating_point_precision","real"); auto& oc = params.sublist("output_control"); oc.set("frequency",1); diff --git a/components/eamxx/src/share/io/tests/output_restart.cpp b/components/eamxx/src/share/io/tests/output_restart.cpp index 9047ed3dfe72..87b584cd9abb 100644 --- a/components/eamxx/src/share/io/tests/output_restart.cpp +++ b/components/eamxx/src/share/io/tests/output_restart.cpp @@ -25,7 +25,7 @@ namespace scream { -constexpr Real FillValue = constants::fill_value; +constexpr Real fill_value = constants::fill_value; std::shared_ptr get_test_fm(const std::shared_ptr& grid); @@ -75,7 +75,6 @@ TEST_CASE("output_restart","io") ekat::ParameterList output_params; output_params.set("floating_point_precision","real"); output_params.set>("field_names",{"field_1", "field_2", "field_3", "field_4","field_5"}); - output_params.set("fill_value",FillValue); output_params.set("flush_frequency",1); output_params.sublist("restart").set("force_new_file",false); output_params.sublist("output_control").set("frequency_units","nsteps"); @@ -272,8 +271,8 @@ void time_advance (const FieldManager& fm, if (fname == "field_5") { // field_5 is used to test restarts w/ filled values, so // we cycle between filled and unfilled states. - v(i,j,k) = (v(i,j,k)==FillValue) ? dt : - ( (v(i,j,k)==1.0) ? 2.0*dt : FillValue ); + v(i,j,k) = (v(i,j,k)==fill_value) ? dt : + ( (v(i,j,k)==1.0) ? 2.0*dt : fill_value ); } else { v(i,j,k) += dt; } diff --git a/components/eamxx/src/share/tests/combine_ops.cpp b/components/eamxx/src/share/tests/combine_ops.cpp index 2dcf49d6ca23..99b6da2003d5 100644 --- a/components/eamxx/src/share/tests/combine_ops.cpp +++ b/components/eamxx/src/share/tests/combine_ops.cpp @@ -20,7 +20,6 @@ TEST_CASE ("combine_ops") { constexpr auto Divide = CombineMode::Divide; constexpr auto Max = CombineMode::Max; constexpr auto Min = CombineMode::Min; - constexpr auto fv_val = constants::fill_value; const pack_type two (2.0); const pack_type four (4.0); @@ -37,7 +36,7 @@ TEST_CASE ("combine_ops") { combine(two,x,2.0,1.0); REQUIRE ( (x==six).all() ); - fill_aware_combine(fv,x,fv_val,2.0,1.0); + combine(fv,x,2.0,1.0); if (not (x==six).all() ) { std::cout << "x: " << x << "\n"; std::cout << " x[0]: " << std::setprecision(18) << x[0] << "\n"; @@ -48,26 +47,26 @@ TEST_CASE ("combine_ops") { x = two; combine(two,x,1,1); REQUIRE ( (x==four).all() ); - fill_aware_combine(fv,x,fv_val,1,1); + combine(fv,x,1,1); REQUIRE ( (x==four).all() ); x = four; combine(two,x,1,1); REQUIRE ( (x==two).all() ); - fill_aware_combine(fv,x,fv_val,1,1); + combine(fv,x,1,1); REQUIRE ( (x==two).all() ); x = two; combine(four,x,1,1); REQUIRE ( (x==four).all() ); - fill_aware_combine(fv,x,fv_val,1,1); + combine(fv,x,1,1); REQUIRE ( (x==four).all() ); x = four; combine(two,x,1,1); REQUIRE ( (x==two).all() ); x = fv_times_ten; - fill_aware_combine(fv,x,fv_val,1,1); + combine(fv,x,1,1); REQUIRE ( (x==fv_times_ten).all() ); } diff --git a/components/eamxx/src/share/tests/field_tests.cpp b/components/eamxx/src/share/tests/field_tests.cpp index 4074b04a4419..3f2ff5308d9b 100644 --- a/components/eamxx/src/share/tests/field_tests.cpp +++ b/components/eamxx/src/share/tests/field_tests.cpp @@ -890,7 +890,7 @@ TEST_CASE ("update") { // Check that updating with rhs==fill_value ignores the rhs f3.deep_copy(constants::fill_value); - f3.get_header().set_extra_data("mask_value",constants::fill_value); + f3.get_header().set_may_be_filled(true); f2.deep_copy(1.0); f2.max(f3); REQUIRE (views_are_equal(f2,one)); @@ -914,7 +914,7 @@ TEST_CASE ("update") { // Check that updating with rhs==fill_value ignores the rhs f3.deep_copy(constants::fill_value); - f3.get_header().set_extra_data("mask_value",constants::fill_value); + f3.get_header().set_may_be_filled(true); f2.deep_copy(1); f2.max(f3); REQUIRE (views_are_equal(f2,one)); @@ -969,7 +969,7 @@ TEST_CASE ("update") { one.deep_copy(1.0); f3.deep_copy(constants::fill_value); - f3.get_header().set_extra_data("mask_value",constants::fill_value); + f3.get_header().set_may_be_filled(true); f2.deep_copy(1.0); f2.update(f3,1,1); if (not views_are_equal(f2,one)) { @@ -1001,7 +1001,7 @@ TEST_CASE ("update") { one.deep_copy(1); f3.deep_copy(constants::fill_value); - f3.get_header().set_extra_data("mask_value",constants::fill_value); + f3.get_header().set_may_be_filled(true); f2.deep_copy(1); f2.update(f3,1,1); REQUIRE (views_are_equal(f2,one)); diff --git a/components/eamxx/src/share/tests/vertical_remapper_tests.cpp b/components/eamxx/src/share/tests/vertical_remapper_tests.cpp index d997587e1ab4..e93e332d9463 100644 --- a/components/eamxx/src/share/tests/vertical_remapper_tests.cpp +++ b/components/eamxx/src/share/tests/vertical_remapper_tests.cpp @@ -15,7 +15,7 @@ constexpr auto Mask = VerticalRemapper::Mask; constexpr auto Top = VerticalRemapper::Top; constexpr auto Bot = VerticalRemapper::Bot; constexpr auto TopBot = VerticalRemapper::TopAndBot; -constexpr Real mask_val = -99999.0; +constexpr Real fill_val = constants::fill_value; template void print (const std::string& fmt, const ekat::Comm& comm, Args&&... args) { @@ -48,7 +48,10 @@ build_grid(const ekat::Comm& comm, const int nldofs, const int nlevs) // Helper function to create fields Field -create_field(const std::string& name, const std::shared_ptr& grid, const bool twod, const bool vec, const bool mid = false, const int ps = 1) +create_field(const std::string& name, + const std::shared_ptr& grid, + const bool create_mask, const bool twod, + const bool vec, const bool mid = false, const int ps = 1) { using namespace ShortFieldTagsNames; constexpr int vec_dim = 3; @@ -62,6 +65,15 @@ create_field(const std::string& name, const std::shared_ptr& Field f(fid); f.get_header().get_alloc_properties().request_allocation(ps); f.allocate_view(); + + if (create_mask) { + // Set a mask (1=filled, 0=valid) to be used in testing the results + FieldIdentifier mask_id("is_filled",fl,units,grid->name(),DataType::IntType); + Field mask(mask_id); + mask.allocate_view(); + mask.deep_copy(0); // Init with "all good" + f.get_header().set_extra_data("is_filled",mask); + } return f; } @@ -161,28 +173,43 @@ void extrapolate (const Field& p_src, const Field& p_tgt, const Field& f, const int ncols = l.dims().front(); const int nlevs = l.dims().back(); const int nlevs_src = p_src.get_header().get_identifier().get_layout().dims().back(); - // print_field_hyperslab(p_src); + auto mask = f.get_header().get_extra_data("is_filled"); + switch (l.type()) { case LayoutType::Scalar2D: break; case LayoutType::Vector2D: break; case LayoutType::Scalar3D: { const auto v = f.get_view(); + const auto m = mask.get_view(); + for (int i=0; ipmax) { - v(i,j) = etype_bot==Mask ? mask_val : data_func(i,0,pmax); + if (etype_bot==Mask) { + v(i,j) = fill_val; + m(i,j) = 1; + } else { + v(i,j) = data_func(i,0,pmax); + } } else if (p(); + const auto m = mask.get_view(); + for (int i=0; ipmax) { - v(i,j,k) = etype_bot==Mask ? mask_val : data_func(i,j,pmax); + if (etype_bot==Mask) { + v(i,j,k) = fill_val; + m(i,j,k) = 1; + } else { + v(i,j,k) = data_func(i,j,pmax); + } } else if (p creating src/tgt pressure fields ... done!\n",comm); print (" -> creating fields ... done!\n",comm); - auto src_s2d = create_field("s2d", src_grid,true,false); - auto src_v2d = create_field("v2d", src_grid,true,true); - auto src_s3d_m = create_field("s3d_m",src_grid,false,false,true, 1); - auto src_s3d_i = create_field("s3d_i",src_grid,false,false,false,SCREAM_PACK_SIZE); - auto src_v3d_m = create_field("v3d_m",src_grid,false,true ,true, 1); - auto src_v3d_i = create_field("v3d_i",src_grid,false,true ,false,SCREAM_PACK_SIZE); - - auto tgt_s2d = create_field("s2d", tgt_grid,true,false); - auto tgt_v2d = create_field("v2d", tgt_grid,true,true); - auto tgt_s3d_m = create_field("s3d_m",tgt_grid,false,false,true, 1); - auto tgt_s3d_i = create_field("s3d_i",tgt_grid,false,false,true, SCREAM_PACK_SIZE); - auto tgt_v3d_m = create_field("v3d_m",tgt_grid,false,true ,true, 1); - auto tgt_v3d_i = create_field("v3d_i",tgt_grid,false,true ,true, SCREAM_PACK_SIZE); - - auto expected_s2d = tgt_s2d.clone(); - auto expected_v2d = tgt_v2d.clone(); - auto expected_s3d_m = tgt_s3d_m.clone(); - auto expected_s3d_i = tgt_s3d_i.clone(); - auto expected_v3d_m = tgt_v3d_m.clone(); - auto expected_v3d_i = tgt_v3d_i.clone(); + auto src_s2d = create_field("s2d", src_grid,false,true,false); + auto src_v2d = create_field("v2d", src_grid,false,true,true); + auto src_s3d_m = create_field("s3d_m",src_grid,false,false,false,true, 1); + auto src_s3d_i = create_field("s3d_i",src_grid,false,false,false,false,SCREAM_PACK_SIZE); + auto src_v3d_m = create_field("v3d_m",src_grid,false,false,true ,true, 1); + auto src_v3d_i = create_field("v3d_i",src_grid,false,false,true ,false,SCREAM_PACK_SIZE); + + auto tgt_s2d = create_field("s2d", tgt_grid,false,true,false); + auto tgt_v2d = create_field("v2d", tgt_grid,false,true,true); + auto tgt_s3d_m = create_field("s3d_m",tgt_grid,false,false,false,true, 1); + auto tgt_s3d_i = create_field("s3d_i",tgt_grid,false,false,false,true, SCREAM_PACK_SIZE); + auto tgt_v3d_m = create_field("v3d_m",tgt_grid,false,false,true ,true, 1); + auto tgt_v3d_i = create_field("v3d_i",tgt_grid,false,false,true ,true, SCREAM_PACK_SIZE); + + // For expected fields, also create the mask extra data + auto expected_s2d = create_field("s2d", tgt_grid,true,true,false); + auto expected_v2d = create_field("v2d", tgt_grid,true,true,true); + auto expected_s3d_m = create_field("s3d_m",tgt_grid,true,false,false,true, 1); + auto expected_s3d_i = create_field("s3d_i",tgt_grid,true,false,false,true, SCREAM_PACK_SIZE); + auto expected_v3d_m = create_field("v3d_m",tgt_grid,true,false,true ,true, 1); + auto expected_v3d_i = create_field("v3d_i",tgt_grid,true,false,true ,true, SCREAM_PACK_SIZE); print (" -> creating fields ... done!\n",comm); // -------------------------------------- // @@ -455,8 +494,6 @@ TEST_CASE ("vertical_remapper") { remap->set_target_pressure (pmid_tgt, pint_tgt); remap->set_extrapolation_type(etype_top,Top); remap->set_extrapolation_type(etype_bot,Bot); - REQUIRE_THROWS (remap->set_mask_value(std::numeric_limits::quiet_NaN())); - remap->set_mask_value(mask_val); // Only needed if top and/or bot use etype=Mask remap->register_field(src_s2d, tgt_s2d); remap->register_field(src_v2d, tgt_v2d); @@ -513,43 +550,25 @@ TEST_CASE ("vertical_remapper") { using namespace Catch::Matchers; Real tol = 10*std::numeric_limits::epsilon(); + auto check = [&](const Field& computed, Field& expected) { + auto diff = computed.clone("diff"); + diff.update(expected,1,-1); + + // Now set expected=0 where it's filled, so we can compute the norm + // (summing x[i]*x[i] would overflow in SP, causing NaN's) + auto mask = expected.get_header().get_extra_data("is_filled"); + expected.deep_copy(0,mask); + auto ex_norm = frobenius_norm(expected); + REQUIRE (frobenius_norm(diff) check tgt fields ...\n",comm); - { - auto diff = tgt_s2d.clone("diff"); - auto ex_norm = frobenius_norm(expected_s2d); - diff.update(expected_s2d,1/ex_norm,-1/ex_norm); - REQUIRE (frobenius_norm(diff)(expected_v2d); - diff.update(expected_v2d,1/ex_norm,-1/ex_norm); - REQUIRE (frobenius_norm(diff)(expected_s3d_m); - diff.update(expected_s3d_m,1/ex_norm,-1/ex_norm); - REQUIRE (frobenius_norm(diff)(expected_s3d_i); - diff.update(expected_s3d_i,1/ex_norm,-1/ex_norm); - REQUIRE (frobenius_norm(diff)(expected_v3d_m); - diff.update(expected_v3d_m,1 / ex_norm,-1 / ex_norm); - REQUIRE (frobenius_norm(diff)(expected_v3d_i); - diff.update(expected_v3d_i,1 / ex_norm,-1 / ex_norm); - REQUIRE (frobenius_norm(diff) check tgt fields ... done!\n",comm); } } diff --git a/components/eamxx/src/share/util/eamxx_combine_ops.hpp b/components/eamxx/src/share/util/eamxx_combine_ops.hpp index f23c48f2f984..e4c5ed478e54 100644 --- a/components/eamxx/src/share/util/eamxx_combine_ops.hpp +++ b/components/eamxx/src/share/util/eamxx_combine_ops.hpp @@ -39,6 +39,8 @@ enum class CombineMode { Min // out = min(beta*out,alpha*in) }; +namespace impl { + // Small helper functions to combine a new value with an old one. // The template argument help reducing the number of operations // performed (the if is resolved at compile time). In the most @@ -46,9 +48,7 @@ enum class CombineMode { // result = (op beta*result alpha*newVal) (where op can be +, *, /, max, min) // This routine should have no overhead compared to a manual // update (assuming you call it with the proper CM) - -template::scalar_type> +template KOKKOS_FORCEINLINE_FUNCTION void combine (const ScalarIn& newVal, ScalarOut& result, const CoeffType alpha, const CoeffType beta) @@ -57,7 +57,7 @@ void combine (const ScalarIn& newVal, ScalarOut& result, using ekat::impl::min; switch (CM) { case CombineMode::Replace: - result = alpha*newVal; + result = newVal; break; case CombineMode::Update: result *= beta; @@ -70,66 +70,52 @@ void combine (const ScalarIn& newVal, ScalarOut& result, result /= (alpha/beta) * newVal; break; case CombineMode::Max: - result = max(beta*result,static_cast(alpha*newVal)); + result = max(beta*result,static_cast(alpha*newVal)); break; case CombineMode::Min: - result = min(beta*result,static_cast(alpha*newVal)); + result = min(beta*result,static_cast(alpha*newVal)); break; + default: + EKAT_KERNEL_ASSERT ("Unsupported/unexpected combine mode.\n"); } } -// 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::scalar_type> -KOKKOS_FORCEINLINE_FUNCTION -std::enable_if_t::is_simd or - ekat::ScalarTraits::is_simd> -fill_aware_combine (const ScalarIn& newVal, ScalarOut& result, - const typename ekat::ScalarTraits::scalar_type fill_val, - const CoeffType alpha, const CoeffType beta) -{ - if constexpr (CM==CombineMode::Replace) { - return combine(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(newVal,tmp,alpha,beta); - where = tmp; - } -} +} // namespace impl -template::scalar_type> KOKKOS_FORCEINLINE_FUNCTION -std::enable_if_t::is_simd and - not ekat::ScalarTraits::is_simd> -fill_aware_combine (const ScalarIn& newVal, ScalarOut& result, - const typename ekat::ScalarTraits::scalar_type fill_val, - const CoeffType alpha, const CoeffType beta) +void combine (const ScalarIn& newVal, ScalarOut& result, + const CoeffType alpha, const CoeffType beta) { - if constexpr (CM==CombineMode::Replace) { - return combine(newVal,result,alpha,beta); + // If not fill-aware, or if CM==Replace, we don't need to check newValue + if constexpr (not fill_aware or CM==CombineMode::Replace) { + return impl::combine(newVal,result,alpha,beta); } - if (newVal!=fill_val) { - combine(newVal,result,alpha,beta); + // For the non-simd type case, we can avoid ekat::where, and simply check newVal against fill_value + if constexpr (ekat::ScalarTraits::is_simd) { + using inner_type = typename ekat::ScalarTraits::scalar_type; + constexpr auto fill_val = constants::fill_value; + auto where = ekat::where(newVal!=fill_val,result); + if (where.any()) { + auto tmp = result; + impl::combine(newVal,tmp,alpha,beta); + where = tmp; + } + } else { + if (newVal!=constants::fill_value) { + impl::combine(newVal,result,alpha,beta); + } } } diff --git a/components/eamxx/src/share/util/eamxx_data_interpolation.cpp b/components/eamxx/src/share/util/eamxx_data_interpolation.cpp index f2eea2743bf7..cb2473631844 100644 --- a/components/eamxx/src/share/util/eamxx_data_interpolation.cpp +++ b/components/eamxx/src/share/util/eamxx_data_interpolation.cpp @@ -539,10 +539,6 @@ create_vert_remapper (const VertRemapData& data) vremap->set_extrapolation_type(s2et(data.extrap_top),VerticalRemapper::Top); vremap->set_extrapolation_type(s2et(data.extrap_bot),VerticalRemapper::Bot); - // Set the mask value only if needed. RemapData has a default that is invalid for VerticalRemapper - if (data.extrap_bot=="Mask" or data.extrap_top=="Mask") { - vremap->set_mask_value(data.mask_value); - } m_vert_remapper = vremap; } diff --git a/components/eamxx/src/share/util/eamxx_data_interpolation.hpp b/components/eamxx/src/share/util/eamxx_data_interpolation.hpp index d292f1b7418b..1511c95906e4 100644 --- a/components/eamxx/src/share/util/eamxx_data_interpolation.hpp +++ b/components/eamxx/src/share/util/eamxx_data_interpolation.hpp @@ -34,7 +34,6 @@ class DataInterpolation VRemapType vr_type = None; std::string extrap_top = "P0"; std::string extrap_bot = "P0"; - Real mask_value = std::numeric_limits::quiet_NaN(); // Unused for P0 extrapolation std::string pname; // What we need to load from nc file Field pmid, pint; // The model pmid/pint std::shared_ptr custom_remapper; // Use this custom remapper diff --git a/components/eamxx/src/share/util/eamxx_time_interpolation.cpp b/components/eamxx/src/share/util/eamxx_time_interpolation.cpp index ba09db3be217..d16efd61440e 100644 --- a/components/eamxx/src/share/util/eamxx_time_interpolation.cpp +++ b/components/eamxx/src/share/util/eamxx_time_interpolation.cpp @@ -164,47 +164,7 @@ void TimeInterpolation::initialize_data_from_files() input_params.set("filename",triplet_curr.filename); m_file_data_atm_input = std::make_shared(input_params,m_fm_time1); m_file_data_atm_input->set_logger(m_logger); - // Assign the mask value gathered from the FillValue found in the source file. - // TODO: Should we make it possible to check if FillValue is in the metadata and only assign mask_value if it is? - for (auto& name : m_field_names) { - auto& field0 = m_fm_time0->get_field(name); - auto& field1 = m_fm_time1->get_field(name); - auto& field_out = m_interp_fields.at(name); - - auto set_fill_value = [&](const auto var_fill_value) { - const auto dt = field_out.data_type(); - if (dt==DataType::FloatType) { - field0.get_header().set_extra_data("mask_value",static_cast(var_fill_value)); - field1.get_header().set_extra_data("mask_value",static_cast(var_fill_value)); - field_out.get_header().set_extra_data("mask_value",static_cast(var_fill_value)); - } else if (dt==DataType::DoubleType) { - field0.get_header().set_extra_data("mask_value",static_cast(var_fill_value)); - field1.get_header().set_extra_data("mask_value",static_cast(var_fill_value)); - field_out.get_header().set_extra_data("mask_value",static_cast(var_fill_value)); - } else { - EKAT_ERROR_MSG ( - "[TimeInterpolation] Unexpected/unsupported field data type.\n" - " - field name: " + field_out.name() + "\n" - " - data type : " + e2str(dt) + "\n"); - } - }; - - const auto& pio_var = scorpio::get_var(triplet_curr.filename,name); - if (scorpio::refine_dtype(pio_var.nc_dtype)=="float") { - auto var_fill_value = scorpio::get_attribute(triplet_curr.filename,name,"_FillValue"); - set_fill_value(var_fill_value); - } else if (scorpio::refine_dtype(pio_var.nc_dtype)=="double") { - auto var_fill_value = scorpio::get_attribute(triplet_curr.filename,name,"_FillValue"); - set_fill_value(var_fill_value); - } else { - EKAT_ERROR_MSG ( - "Unrecognized/unsupported data type\n" - " - filename: " + triplet_curr.filename + "\n" - " - varname : " + name + "\n" - " - dtype : " + pio_var.dtype + "\n"); - } - } // Read first snap of data and shift to time0 read_data(); shift_data(); @@ -351,24 +311,6 @@ void TimeInterpolation::read_data() input_params.set("filename",triplet_curr.filename); m_file_data_atm_input = std::make_shared(input_params,m_fm_time1); m_file_data_atm_input->set_logger(m_logger); - // Also determine the FillValue, if used - // TODO: Should we make it possible to check if FillValue is in the metadata and only assign mask_value if it is? - for (auto& name : m_field_names) { - auto& field = m_fm_time1->get_field(name); - const auto dt = field.data_type(); - if (dt==DataType::FloatType) { - auto var_fill_value = scorpio::get_attribute(triplet_curr.filename,name,"_FillValue"); - field.get_header().set_extra_data("mask_value",var_fill_value); - } else if (dt==DataType::DoubleType) { - auto var_fill_value = scorpio::get_attribute(triplet_curr.filename,name,"_FillValue"); - field.get_header().set_extra_data("mask_value",var_fill_value); - } else { - EKAT_ERROR_MSG ( - "[TimeInterpolation] Unexpected/unsupported field data type.\n" - " - field name: " + field.name() + "\n" - " - data type : " + e2str(dt) + "\n"); - } - } } if (m_logger) { diff --git a/components/eamxx/tests/single-process/surface_coupling/output.yaml b/components/eamxx/tests/single-process/surface_coupling/output.yaml index 837d5bb924cc..7b0613ed050a 100644 --- a/components/eamxx/tests/single-process/surface_coupling/output.yaml +++ b/components/eamxx/tests/single-process/surface_coupling/output.yaml @@ -2,7 +2,6 @@ --- filename_prefix: surface_coupling_output averaging_type: instant -fill_value: -99999.0 field_names: # IMPORTER - sfc_alb_dir_vis