From be93f5d2ec7a46e55ea7d2fa21d59eecdd62e7bd Mon Sep 17 00:00:00 2001 From: "Michael J. Schmidt" Date: Thu, 21 Mar 2024 17:04:14 -0600 Subject: [PATCH 01/16] multi-slice subfield written and compiling --- components/eamxx/src/share/field/field.cpp | 68 ++++++++++++++++++- components/eamxx/src/share/field/field.hpp | 15 +++- .../src/share/field/field_alloc_prop.cpp | 52 ++++++++++++++ .../src/share/field/field_alloc_prop.hpp | 44 ++++++++---- .../eamxx/src/share/field/field_header.cpp | 37 ++++++++-- .../eamxx/src/share/field/field_header.hpp | 11 +++ 6 files changed, 206 insertions(+), 21 deletions(-) diff --git a/components/eamxx/src/share/field/field.cpp b/components/eamxx/src/share/field/field.cpp index e022278d503..11c8e695b4b 100644 --- a/components/eamxx/src/share/field/field.cpp +++ b/components/eamxx/src/share/field/field.cpp @@ -82,6 +82,7 @@ Field Field:: subfield (const std::string& sf_name, const ekat::units::Units& sf_units, const int idim, const int index, const bool dynamic) const { + // const auto& id = m_header->get_identifier(); const auto& lt = id.get_layout(); @@ -115,6 +116,50 @@ subfield (const int idim, const int index, const bool dynamic) const { return subfield(m_header->get_identifier().name(),idim,index,dynamic); } +// slice at index idim, entries \in [index_beg, index_end] +Field Field::subfield(const std::string& sf_name, + const ekat::units::Units& sf_units, const int idim, + const int index_beg, const int index_end, + const bool dynamic) const { + + const auto& id = m_header->get_identifier(); + const auto& lt = id.get_layout(); + + // Sanity checks + EKAT_REQUIRE_MSG( + is_allocated(), + "Error! Input field must be allocated in order to subview it.\n"); + EKAT_REQUIRE_MSG(idim == 0 || idim == 1, + "Error! Subview dimension index must be either 0 or 1.\n"); + + auto sf_layout = + lt.clone_with_different_extent(idim, index_end - index_beg + 1); + // Create identifier for subfield + FieldIdentifier sf_id(sf_name, sf_layout, sf_units, id.get_grid_name()); + + // Create empty subfield, then set header and views + // Note: we can access protected members, since it's the same type + Field sf; + sf.m_header = create_subfield_header(sf_id, m_header, idim, index_beg, + index_end, dynamic); + sf.m_data = m_data; + + return sf; +} + +Field Field::subfield(const std::string& sf_name, const int idim, + const int index_beg, const int index_end, + const bool dynamic) const { + const auto& id = m_header->get_identifier(); + return subfield(sf_name, id.get_units(), idim, index_beg, index_end, dynamic); +} + +Field Field::subfield(const int idim, const int index_beg, const int index_end, + const bool dynamic) const { + return subfield(m_header->get_identifier().name(), idim, index_beg, index_end, + dynamic); +} + Field Field:: get_component (const int i, const bool dynamic) { const auto& layout = get_header().get_identifier().get_layout(); @@ -127,11 +172,32 @@ get_component (const int i, const bool dynamic) { EKAT_REQUIRE_MSG (i>=0 && i= 0 && i2 < layout.dim(idim), + "Error! Component index range out of bounds [0," + + std::to_string(layout.dim(idim)) + ").\n"); + EKAT_REQUIRE_MSG(i1 < i2, "Error! Invalid component indices (i1 >= i2).\n"); + + // Add _$i1-$i2 to the field name, to avoid issues if the subfield is stored + // in some structure that requires unique names (e.g., a remapper) + return subfield(fname + "_" + std::to_string(i1) + "-" + std::to_string(i2), + idim, i1, i2, dynamic); +} + bool Field::equivalent(const Field& rhs) const { return (m_header==rhs.m_header && diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index 4a4b20525a3..c5c7ef6de40 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -241,11 +241,24 @@ class Field { Field subfield (const std::string& sf_name, const int idim, const int index, const bool dynamic = false) const; Field subfield (const int idim, const int k, const bool dynamic = false) const; - + // get subfield to extract multiple slices in a continuous range of indices + // e.g., (in matlab syntax) subf = f.subfield(:, 1:3, :) + // but NOT subf = f.subfield(:, [1, 3, 4], :) + // NOTE: use std::pair(index_beg, index_end) or a kokkos::pair()? + Field subfield (const std::string& sf_name, const ekat::units::Units& sf_units, + const int idim, const int index_beg, const int index_end, + const bool dynamic = false) const; + Field subfield (const std::string& sf_name, const int idim, + const int index_beg, const int index_end, + const bool dynamic = false) const; + Field subfield (const int idim, const int k_beg, const int k_end, + const bool dynamic = false) const; // If this field is a vector field, get a subfield for the ith component. // If dynamic = true, it is possible to "reset" the component index at runtime. // Note: throws if this is not a vector field. Field get_component (const int i, const bool dynamic = false); + // version for slicing across multiple, contiguous indices + Field get_component (const int i1, const int i2, const bool dynamic = false); // Checks whether the underlying view has been already allocated. bool is_allocated () const { return m_data.d_view.data()!=nullptr; } diff --git a/components/eamxx/src/share/field/field_alloc_prop.cpp b/components/eamxx/src/share/field/field_alloc_prop.cpp index fcffc821bba..f2375869f84 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.cpp +++ b/components/eamxx/src/share/field/field_alloc_prop.cpp @@ -77,6 +77,58 @@ subview (const int idim, const int k, const bool dynamic) const { return props; } +FieldAllocProp FieldAllocProp::subview(const int idim, const int k_beg, + const int k_end, + const bool dynamic) const { + EKAT_REQUIRE_MSG( + is_committed(), + "Error! Subview requires alloc properties to be committed.\n"); + EKAT_REQUIRE_MSG( + idim == 0 || idim == 1, + "Error! Subviewing is only allowed along first or second dimension.\n"); + EKAT_REQUIRE_MSG(idim < m_layout.rank(), + "Error! Dimension index out of bounds.\n"); + EKAT_REQUIRE_MSG(k_beg <= k_end, + "Error! Slice indices are invalid (non-increasing).\n"); + EKAT_REQUIRE_MSG( + k_beg >= 0 && k_end < m_layout.dim(idim), + "Error! Index range along the dimension is out of bounds.\n"); + + // Set new layout basic stuff + FieldAllocProp props(m_scalar_type_size); + props.m_committed = true; + props.m_scalar_type_size = m_scalar_type_size; + props.m_layout = + m_layout.clone_with_different_extent(idim, k_end - k_beg + 1); + + // Output is contiguous if either + // - this->m_contiguous=true AND idim==0 + // - m_layout.dim(i)==1 for all i. * Let's say you have one class that uses the field as View, - * and another class that used the field in a packed way, taht is, - * as View**>. The second view will have dimensions (3,4), + * and another class that used the field in a packed way, that is, + * as View**>. The second view will have dimensions (3,3), in terms of its value type (i.e., Pack). * This means the field needs an allocation bigger than 30 Real's, namely - * 3*4*4=36 Real's. This class does the book-keeping for the allocation size, + * 3*4*3=36 Real's. This class does the book-keeping for the allocation size, * so that 1) the field can be allocated with enough memory to accommodate * all requests, 2) customers of the field can check what the allocation * is, so that they know whether there is padding in the field, and @@ -61,16 +61,30 @@ namespace scream // Helper struct to store info related to the subviewing process struct SubviewInfo { - SubviewInfo () = default; - SubviewInfo (const int dim, const int slice, const int extent, const bool is_dynamic) : - dim_idx(dim), slice_idx(slice), dim_extent(extent), dynamic(is_dynamic) {} - SubviewInfo (const SubviewInfo&) = default; - SubviewInfo& operator= (const SubviewInfo&) = default; - - int dim_idx = -1; // Dimension along which slicing happened - int slice_idx = -1; // Slice along dimension $dim_idx - int dim_extent = -1; // Extent of dimension $dim_idx - bool dynamic = false; // Whether this is a dynamic subview (slice_idx can change) + SubviewInfo() = default; + SubviewInfo(const int dim, const int slice, const int extent, + const bool is_dynamic) + : dim_idx(dim), slice_idx(slice), dim_extent(extent), + dynamic(is_dynamic) {} + // multi-slice subview across contiguous indices + SubviewInfo(const int dim, const int slice_beg, const int slice_end, + const int extent, const bool is_dynamic) + : dim_idx(dim), slice_idx_beg(slice_beg), slice_idx_end(slice_end), + dim_extent(extent), dynamic(is_dynamic) {} + SubviewInfo(const SubviewInfo&) = default; + SubviewInfo& operator=(const SubviewInfo&) = default; + + // TODO: is it better to do it this way, or to create an alternative struct? + int dim_idx = -1; // Dimension along which slicing happened + // Slice along dimension $dim_idx (uninitialized if multi-slice subview) + int slice_idx = -1; + // beginning slice index for multi-slice (remains == -1 if single-slice subview) + int slice_idx_beg = -1; + // ending slice index for multi-slice (remains == -1 if single-slice subview) + int slice_idx_end = -1; + int dim_extent = -1; // Extent of dimension $dim_idx + bool dynamic = + false; // Whether this is a dynamic subview (slice_idx can change) }; inline bool operator== (const SubviewInfo& lhs, const SubviewInfo& rhs) { @@ -95,6 +109,10 @@ class FieldAllocProp { // possible to change the entry index (k) at runtime. FieldAllocProp subview (const int idim, const int k, const bool dynamic) const; + // multi-slice subview over contiguous indices + FieldAllocProp subview (const int idim, const int k_beg, const int k_end, + const bool dynamic) const; + // Request allocation able to accommodate a pack of ScalarType of the given pack size void request_allocation (const int pack_size = 1); diff --git a/components/eamxx/src/share/field/field_header.cpp b/components/eamxx/src/share/field/field_header.cpp index 9bb642f2236..934084b3080 100644 --- a/components/eamxx/src/share/field/field_header.cpp +++ b/components/eamxx/src/share/field/field_header.cpp @@ -2,8 +2,7 @@ #include "ekat/std_meta/ekat_std_utils.hpp" -namespace scream -{ +namespace scream { FieldHeader::FieldHeader (const identifier_type& id) : m_identifier (id) @@ -33,10 +32,7 @@ set_extra_data (const std::string& key, } } -std::shared_ptr -FieldHeader:: -alias (const std::string& name) const -{ +std::shared_ptr FieldHeader::alias(const std::string& name) const { auto fh = create_header(get_identifier().alias(name)); fh->m_tracking = m_tracking; fh->m_alloc_prop = m_alloc_prop; @@ -72,4 +68,33 @@ create_subfield_header (const FieldIdentifier& id, return fh; } +// subfield with multiple, contiguous slices +std::shared_ptr +create_subfield_header(const FieldIdentifier& id, + std::shared_ptr parent, const int idim, + const int k_beg, const int k_end, const bool dynamic) { + // Sanity checks + EKAT_REQUIRE_MSG(parent != nullptr, + "Error! Invalid pointer for parent header.\n"); + EKAT_REQUIRE_MSG(k_end > k_beg, + "Error! Slice indices are invalid (non-increasing).\n"); + + // Create header, and set up parent/child + auto fh = create_header(id); + fh->create_parent_child_link(parent); + + // Create tracking, and set up parent/child + fh->m_tracking = create_tracking(); + fh->m_tracking->create_parent_child_link(parent->get_tracking_ptr()); + if (parent->get_tracking().get_time_stamp().is_valid()) { + fh->m_tracking->update_time_stamp(parent->get_tracking().get_time_stamp()); + } + + // Create alloc props + fh->m_alloc_prop = std::make_shared( + parent->get_alloc_properties().subview(idim, k_beg, k_end, dynamic)); + + return fh; +} + } // namespace scream diff --git a/components/eamxx/src/share/field/field_header.hpp b/components/eamxx/src/share/field/field_header.hpp index 511c6b6f505..84284272013 100644 --- a/components/eamxx/src/share/field/field_header.hpp +++ b/components/eamxx/src/share/field/field_header.hpp @@ -91,6 +91,12 @@ class FieldHeader : public FamilyTracking { create_subfield_header (const FieldIdentifier&, std::shared_ptr, const int, const int, const bool); + // for creating multi-slice subfield (continuous indices) + friend std::shared_ptr + create_subfield_header (const FieldIdentifier&, + std::shared_ptr, + const int idim, const int k_beg, const int k_end, + const bool dynamic); // NOTE: the identifier *cannot* be a shared_ptr, b/c we // don't foresee sharing an identifier between two @@ -169,6 +175,11 @@ std::shared_ptr create_subfield_header (const FieldIdentifier& id, std::shared_ptr parent, const int idim, const int k, const bool dynamic); +std::shared_ptr +create_subfield_header (const FieldIdentifier& id, + std::shared_ptr parent, + const int idim, const int k_beg, const int k_end, + const bool dynamic); } // namespace scream From ccc846a30d87ddb3caf066d4134145e69dcdbcb4 Mon Sep 17 00:00:00 2001 From: "Michael J. Schmidt" Date: Mon, 1 Apr 2024 17:26:19 -0600 Subject: [PATCH 02/16] got get_strided_view() working for multi-sliced field--ugly impl at the moment --- components/eamxx/src/share/field/field.hpp | 6 +- .../src/share/field/field_alloc_prop.hpp | 8 +- .../eamxx/src/share/field/field_impl.hpp | 89 ++++++++++++++++++- .../eamxx/src/share/tests/field_tests.cpp | 49 ++++++++++ 4 files changed, 144 insertions(+), 8 deletions(-) diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index c5c7ef6de40..b716572f0f4 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -136,6 +136,10 @@ class Field { get_strided_view_type get_strided_view () const; + template + Kokkos::View + get_strided_view (bool special) const; + // These two getters are convenience function for commonly accessed metadata. // The same info can be extracted from the metadata stored in the FieldHeader DataType data_type () const { return get_header().get_identifier().data_type(); } @@ -327,7 +331,7 @@ class Field { auto get_ND_view () const -> if_t<(N,HD>>; - // Metadata (name, rank, dims, customere/providers, time stamp, ...) + // Metadata (name, rank, dims, customer/providers, time stamp, ...) std::shared_ptr m_header; // Actual data. diff --git a/components/eamxx/src/share/field/field_alloc_prop.hpp b/components/eamxx/src/share/field/field_alloc_prop.hpp index c3f870b4280..6f23d579d89 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.hpp +++ b/components/eamxx/src/share/field/field_alloc_prop.hpp @@ -69,14 +69,15 @@ struct SubviewInfo { // multi-slice subview across contiguous indices SubviewInfo(const int dim, const int slice_beg, const int slice_end, const int extent, const bool is_dynamic) - : dim_idx(dim), slice_idx_beg(slice_beg), slice_idx_end(slice_end), + : dim_idx(dim), slice_idx(slice_beg), slice_idx_end(slice_end), dim_extent(extent), dynamic(is_dynamic) {} SubviewInfo(const SubviewInfo&) = default; SubviewInfo& operator=(const SubviewInfo&) = default; - // TODO: is it better to do it this way, or to create an alternative struct? int dim_idx = -1; // Dimension along which slicing happened - // Slice along dimension $dim_idx (uninitialized if multi-slice subview) + // Slice along dimension $dim_idx if taking single-index slice for subfield + // If taking multi-slice subfield, then this is the starting slice index + // e.g., slicing (:, slice_idx_end : slice_idx_end, :) int slice_idx = -1; // beginning slice index for multi-slice (remains == -1 if single-slice subview) int slice_idx_beg = -1; @@ -90,6 +91,7 @@ struct SubviewInfo { inline bool operator== (const SubviewInfo& lhs, const SubviewInfo& rhs) { return lhs.dim_idx==rhs.dim_idx && lhs.slice_idx==rhs.slice_idx && + // FIXME: change for multi-slice (throw error to check?) lhs.dim_extent==rhs.dim_extent && lhs.dynamic==rhs.dynamic; } diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index f192d70bb68..7ad4e231bd8 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -40,7 +40,7 @@ Field (const identifier_type& id, " - field name: " + id.name() + "\n" " - idim: " + std::to_string(i) + "\n" " - layout i-th dim: " + std::to_string(fl.dims()[i]) + "\n" - " - view i-th dim: " + std::to_string(view_d.extent(i)) + "\n"); + " - view i-th dim: " + std::to_string(view_d.extent(i)) + "\n"); } auto& alloc_prop = m_header->get_alloc_properties(); @@ -49,7 +49,7 @@ Field (const identifier_type& id, "Error! Input view has the wrong last extent.\n" " - field name: " + id.name() + "\n" " - layout last dim: " + std::to_string(fl.dims()[N-1]) + "\n" - " - view last dim: " + std::to_string(view_d.extent(N-1)) + "\n"); + " - view last dim: " + std::to_string(view_d.extent(N-1)) + "\n"); // We have a padded view. We don't know what the pack size was, so we pick the largest // power of 2 that divides the last extent @@ -73,8 +73,6 @@ Field (const identifier_type& id, // Create an unmanaged dev view, and its host mirror const auto view_dim = alloc_prop.get_alloc_size(); char* data = reinterpret_cast(view_d.data()); - std::cout << "fl: " << fl.to_string() << "\n" - << "view dim: " << view_dim << "\n"; m_data.d_view = decltype(m_data.d_view)(data,view_dim); m_data.h_view = Kokkos::create_mirror_view(m_data.d_view); @@ -101,6 +99,10 @@ auto Field::get_view () const // Make sure input field is allocated EKAT_REQUIRE_MSG(is_allocated(), "Error! Cannot extract a field's view before allocation happens.\n"); + // FIXME: add check + // EKAT_REQUIRE_MSG(true /*is_multiSlice_subview()*/, + // "Error! Multi-sliced subfield is incompatible--must employ " + // "get_strided_view().\n") EKAT_REQUIRE_MSG (not m_is_read_only || std::is_const::value, "Error! Cannot get a view to non-const data if the field is read-only.\n"); @@ -205,6 +207,83 @@ auto Field::get_strided_view () const return DstView(get_ND_view()); } +template +Kokkos::View Field::get_strided_view (bool special) const +{ + // The destination view type on correct mem space + using DstView = get_strided_view_type; + // The dst value types + using DstValueType = typename DstView::traits::value_type; + // We only allow to reshape to a view of the correct rank + constexpr int DstRank = DstView::rank; + constexpr int DstRankDynamic = DstView::rank_dynamic; + + // Get src details + const auto& alloc_prop = m_header->get_alloc_properties(); + const auto& fl = m_header->get_identifier().get_layout(); + + // Checks + // EKAT_REQUIRE_MSG (DstRank==1 && fl.rank()==1, + // "Error! Strided view only available for rank-1 fields.\n"); + // EKAT_REQUIRE_MSG (DstRankDynamic==1, + // "Error! Strided view not allowed with compile-time dimensions.\n"); + EKAT_REQUIRE_MSG(is_allocated(), + "Error! Cannot extract a field's view before allocation happens.\n"); + EKAT_REQUIRE_MSG (not m_is_read_only || std::is_const::value, + "Error! Cannot get a view to non-const data if the field is read-only.\n"); + EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), + "Error! Source field allocation is not compatible with the requested value type.\n"); + + // Check if this field is a subview of another field + const auto parent = m_header->get_parent().lock(); + if (parent!=nullptr) { + // Parent field has correct layout to reinterpret the view into N+1-dim view + // So create the parent field on the fly, use it to get the N+1-dim view, then subview it. + // NOTE: we can set protected members, since f is the same type of this class. + Field f; + f.m_header = parent; + f.m_data = m_data; + + // // Take 2 dimensional view with normal LayoutRight + // auto v_np1 = f.get_ND_view(); + // const auto vfs_rank = f.get_header().get_identifier().get_layout().rank(); + auto v_fullsize = f.get_ND_view(); + + // Now we can subview v_np1 at the correct slice + const auto& info = m_header->get_alloc_properties().get_subview_info(); + const int idim = info.dim_idx; + const int k = info.slice_idx; + const int k_end = info.slice_idx_end; + + // // So far we can only subview at first or second dimension. + // EKAT_REQUIRE_MSG (idim==0 || idim==1, + // "Error! Subview dimension index is out of bounds.\n"); + + // Use correct subview utility + if (idim==0) { + // FIXME: what's a better way to do this? can't use v_fullsize after the + // logic block b/c of scoping, and it's also not easy to know what type of + // exotic + // we know it's a subview, so if it has the same rank as its parent, then + // we know that it's a multi-slice subview + if (fl.rank() == f.get_header().get_identifier().get_layout().rank()) { + auto svs = ekat::subview(v_fullsize, Kokkos::make_pair(k, k_end), idim); + // for (size_t i = 0; i < 4; i++) + // { + // std::cout << "i = " << i << "\n"; + // std::cout << "svs.stride(i) = " << svs.stride(i) << "\n"; + // std::cout << "svs.extent(i) = " << svs.extent(i) << "\n"; + // } + return svs; + } else { + // return DstView(ekat::subview(v_np1,k)); + } + } else { + // return DstView(ekat::subview_1(v_np1,k)); + } + } +} + template void Field:: deep_copy (const Field& src) { @@ -764,6 +843,7 @@ auto Field::get_ND_view () const -> return ret_type (ptr,kl); } +// TODO: will need to set this up for multi-sliced subfield template auto Field::get_ND_view () const -> if_t,HD>> @@ -773,6 +853,7 @@ auto Field::get_ND_view () const -> "Error! Input Rank must either be 1 (flat array) or the actual field rank.\n"); // Given that N==MaxRank, this field cannot be a subview of another field + // NOTE: this will not be true for multi-slice EKAT_REQUIRE_MSG (m_header->get_parent().expired(), "Error! A view of rank " + std::to_string(MaxRank) + " should not be the subview of another field.\n"); diff --git a/components/eamxx/src/share/tests/field_tests.cpp b/components/eamxx/src/share/tests/field_tests.cpp index a1336ae346f..0411755dda4 100644 --- a/components/eamxx/src/share/tests/field_tests.cpp +++ b/components/eamxx/src/share/tests/field_tests.cpp @@ -308,6 +308,55 @@ TEST_CASE("field", "") { REQUIRE (f3.is_read_only()); } + // Subfields (multi-sliced) + SECTION ("subfield--multi-slice") { + std::vector t1 = {COL,CMP,CMP,LEV}; + std::vector d1 = {5,10,2,24}; + + FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); + + Field f1(fid1); + f1.allocate_view(); + randomize(f1,engine,pdf); + + const int idim = 0; + const int ivar = 2; + const int sl_beg = 2; + const int sl_end = 4; + + auto f3 = f1.subfield(idim, sl_beg, sl_end); + + // const auto fl3 = f3.get_header().get_identifier().get_layout(); + // const auto alprop3 = f3.get_header().get_alloc_properties(); + // std::cout << "fl3.rank() = " << fl3.rank() << "\n"; + // auto f3dim = fl3.dims(); + // for (size_t i = 0; i < 4; i++) + // { + // std::cout << "f3dim(i) = " << f3dim[i] << "\n"; + // } + // auto f3ext = fl3.extents(); + // for (size_t i = 0; i < 4; i++) + // { + // std::cout << "f3ext(i) = " << f3ext(i) << "\n"; + // } + // std::cout << "alprop3.is_subfield() = " << alprop3.is_subfield() << "\n"; + + auto v3_h = f3.get_strided_view(true); + auto v4d_h = f1.get_view(); + + for (size_t i = sl_beg; i < sl_end; i++) { + for (size_t j = 0; j < d1[1]; j++) { + for (size_t k = 0; k < d1[2]; k++) { + for (size_t l = 0; l < d1[3]; l++) { + // std::cout << "v4d_h(i, j, k, l) = " << v4d_h(i, j, k, l) << "\n"; + // std::cout << "v3_h(i - sl_beg, j, k, l) = " << v3_h(i - sl_beg, j, k, l) << "\n"; + REQUIRE(v4d_h(i, j, k, l) == v3_h(i - sl_beg, j, k, l)); + } + } + } + } + } + // Dynamic Subfields SECTION ("dynamic_subfield") { const int vec_dim = 10; From da94fdb73082dce84474ea6a55b74f8c9cf25f4c Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Wed, 3 Apr 2024 19:40:44 -0600 Subject: [PATCH 03/16] multi-sliced subfield/subviews working --- components/eamxx/src/share/field/field.cpp | 4 +- components/eamxx/src/share/field/field.hpp | 11 +-- .../src/share/field/field_alloc_prop.cpp | 12 ++-- .../src/share/field/field_alloc_prop.hpp | 2 +- .../eamxx/src/share/field/field_impl.hpp | 72 +++++-------------- .../eamxx/src/share/tests/field_tests.cpp | 53 ++++++-------- 6 files changed, 56 insertions(+), 98 deletions(-) diff --git a/components/eamxx/src/share/field/field.cpp b/components/eamxx/src/share/field/field.cpp index 11c8e695b4b..1e790ae3cea 100644 --- a/components/eamxx/src/share/field/field.cpp +++ b/components/eamxx/src/share/field/field.cpp @@ -129,8 +129,8 @@ Field Field::subfield(const std::string& sf_name, EKAT_REQUIRE_MSG( is_allocated(), "Error! Input field must be allocated in order to subview it.\n"); - EKAT_REQUIRE_MSG(idim == 0 || idim == 1, - "Error! Subview dimension index must be either 0 or 1.\n"); + // EKAT_REQUIRE_MSG(idim == 0 || idim == 1, + // "Error! Subview dimension index must be either 0 or 1.\n"); auto sf_layout = lt.clone_with_different_extent(idim, index_end - index_beg + 1); diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index b716572f0f4..98169ece28d 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -136,9 +136,11 @@ class Field { get_strided_view_type get_strided_view () const; - template - Kokkos::View - get_strided_view (bool special) const; + // this is the same as above but for a multi-sliced view + // which is to say, also a strided view + template + auto get_multi_sliced_view () const + -> get_strided_view_type, HD>; // These two getters are convenience function for commonly accessed metadata. // The same info can be extracted from the metadata stored in the FieldHeader @@ -245,10 +247,9 @@ class Field { Field subfield (const std::string& sf_name, const int idim, const int index, const bool dynamic = false) const; Field subfield (const int idim, const int k, const bool dynamic = false) const; - // get subfield to extract multiple slices in a continuous range of indices + // subfield fxn to extract multiple slices in a continuous range of indices // e.g., (in matlab syntax) subf = f.subfield(:, 1:3, :) // but NOT subf = f.subfield(:, [1, 3, 4], :) - // NOTE: use std::pair(index_beg, index_end) or a kokkos::pair()? Field subfield (const std::string& sf_name, const ekat::units::Units& sf_units, const int idim, const int index_beg, const int index_end, const bool dynamic = false) const; diff --git a/components/eamxx/src/share/field/field_alloc_prop.cpp b/components/eamxx/src/share/field/field_alloc_prop.cpp index f2375869f84..08d850a9325 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.cpp +++ b/components/eamxx/src/share/field/field_alloc_prop.cpp @@ -77,22 +77,20 @@ subview (const int idim, const int k, const bool dynamic) const { return props; } -FieldAllocProp FieldAllocProp::subview(const int idim, const int k_beg, +FieldAllocProp FieldAllocProp::subview(const int idim, + const int k_beg, const int k_end, const bool dynamic) const { EKAT_REQUIRE_MSG( is_committed(), "Error! Subview requires alloc properties to be committed.\n"); - EKAT_REQUIRE_MSG( - idim == 0 || idim == 1, - "Error! Subviewing is only allowed along first or second dimension.\n"); - EKAT_REQUIRE_MSG(idim < m_layout.rank(), + EKAT_REQUIRE_MSG(idim <= m_layout.rank(), "Error! Dimension index out of bounds.\n"); - EKAT_REQUIRE_MSG(k_beg <= k_end, + EKAT_REQUIRE_MSG(k_beg < k_end, "Error! Slice indices are invalid (non-increasing).\n"); EKAT_REQUIRE_MSG( k_beg >= 0 && k_end < m_layout.dim(idim), - "Error! Index range along the dimension is out of bounds.\n"); + "Error! Slice index range along the idim dimension is out of bounds.\n"); // Set new layout basic stuff FieldAllocProp props(m_scalar_type_size); diff --git a/components/eamxx/src/share/field/field_alloc_prop.hpp b/components/eamxx/src/share/field/field_alloc_prop.hpp index 6f23d579d89..481182fbdc6 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.hpp +++ b/components/eamxx/src/share/field/field_alloc_prop.hpp @@ -106,7 +106,7 @@ class FieldAllocProp { FieldAllocProp& operator= (const FieldAllocProp&); - // Return allocation props of an unmanaged subivew of this field + // Return allocation props of an unmanaged subview of this field // at entry k along dimension idim. If dynamic = true, then it will be // possible to change the entry index (k) at runtime. FieldAllocProp subview (const int idim, const int k, const bool dynamic) const; diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index 7ad4e231bd8..b938cb6070f 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -99,10 +99,6 @@ auto Field::get_view () const // Make sure input field is allocated EKAT_REQUIRE_MSG(is_allocated(), "Error! Cannot extract a field's view before allocation happens.\n"); - // FIXME: add check - // EKAT_REQUIRE_MSG(true /*is_multiSlice_subview()*/, - // "Error! Multi-sliced subfield is incompatible--must employ " - // "get_strided_view().\n") EKAT_REQUIRE_MSG (not m_is_read_only || std::is_const::value, "Error! Cannot get a view to non-const data if the field is read-only.\n"); @@ -207,26 +203,26 @@ auto Field::get_strided_view () const return DstView(get_ND_view()); } -template -Kokkos::View Field::get_strided_view (bool special) const +// NOTE: multi-slicing a view is only supported for strided view return type +template +auto Field::get_multi_sliced_view () const + -> get_strided_view_type, HD> { // The destination view type on correct mem space - using DstView = get_strided_view_type; + using DstView = get_strided_view_type, HD>; // The dst value types using DstValueType = typename DstView::traits::value_type; - // We only allow to reshape to a view of the correct rank - constexpr int DstRank = DstView::rank; - constexpr int DstRankDynamic = DstView::rank_dynamic; // Get src details const auto& alloc_prop = m_header->get_alloc_properties(); const auto& fl = m_header->get_identifier().get_layout(); // Checks - // EKAT_REQUIRE_MSG (DstRank==1 && fl.rank()==1, - // "Error! Strided view only available for rank-1 fields.\n"); - // EKAT_REQUIRE_MSG (DstRankDynamic==1, + // TODO: decide whether a dynamic-rank view is ok + // EKAT_REQUIRE_MSG (DstRankDynamic == 1, // "Error! Strided view not allowed with compile-time dimensions.\n"); + EKAT_REQUIRE_MSG (N == fl.rank(), + "Error! Input Rank must be equal to parent view's rank for multi-sliced subview.\n"); EKAT_REQUIRE_MSG(is_allocated(), "Error! Cannot extract a field's view before allocation happens.\n"); EKAT_REQUIRE_MSG (not m_is_read_only || std::is_const::value, @@ -236,51 +232,23 @@ Kokkos::View Field::get_strided_view (bool speci // Check if this field is a subview of another field const auto parent = m_header->get_parent().lock(); - if (parent!=nullptr) { - // Parent field has correct layout to reinterpret the view into N+1-dim view - // So create the parent field on the fly, use it to get the N+1-dim view, then subview it. - // NOTE: we can set protected members, since f is the same type of this class. + if (parent != nullptr) { Field f; f.m_header = parent; f.m_data = m_data; - // // Take 2 dimensional view with normal LayoutRight - // auto v_np1 = f.get_ND_view(); - // const auto vfs_rank = f.get_header().get_identifier().get_layout().rank(); - auto v_fullsize = f.get_ND_view(); + auto v_fullsize = f.get_ND_view(); // Now we can subview v_np1 at the correct slice const auto& info = m_header->get_alloc_properties().get_subview_info(); - const int idim = info.dim_idx; - const int k = info.slice_idx; - const int k_end = info.slice_idx_end; - - // // So far we can only subview at first or second dimension. - // EKAT_REQUIRE_MSG (idim==0 || idim==1, - // "Error! Subview dimension index is out of bounds.\n"); - - // Use correct subview utility - if (idim==0) { - // FIXME: what's a better way to do this? can't use v_fullsize after the - // logic block b/c of scoping, and it's also not easy to know what type of - // exotic - // we know it's a subview, so if it has the same rank as its parent, then - // we know that it's a multi-slice subview - if (fl.rank() == f.get_header().get_identifier().get_layout().rank()) { - auto svs = ekat::subview(v_fullsize, Kokkos::make_pair(k, k_end), idim); - // for (size_t i = 0; i < 4; i++) - // { - // std::cout << "i = " << i << "\n"; - // std::cout << "svs.stride(i) = " << svs.stride(i) << "\n"; - // std::cout << "svs.extent(i) = " << svs.extent(i) << "\n"; - // } - return svs; - } else { - // return DstView(ekat::subview(v_np1,k)); - } - } else { - // return DstView(ekat::subview_1(v_np1,k)); - } + const int idim = info.dim_idx; + const int k = info.slice_idx; + const int k_end = info.slice_idx_end; + + // this version of ekat::subview overloaded conveniently, so only the single + // version of the call is required here + return DstView(ekat::subview(v_fullsize, + Kokkos::make_pair(k, k_end), idim)); } } @@ -843,7 +811,6 @@ auto Field::get_ND_view () const -> return ret_type (ptr,kl); } -// TODO: will need to set this up for multi-sliced subfield template auto Field::get_ND_view () const -> if_t,HD>> @@ -853,7 +820,6 @@ auto Field::get_ND_view () const -> "Error! Input Rank must either be 1 (flat array) or the actual field rank.\n"); // Given that N==MaxRank, this field cannot be a subview of another field - // NOTE: this will not be true for multi-slice EKAT_REQUIRE_MSG (m_header->get_parent().expired(), "Error! A view of rank " + std::to_string(MaxRank) + " should not be the subview of another field.\n"); diff --git a/components/eamxx/src/share/tests/field_tests.cpp b/components/eamxx/src/share/tests/field_tests.cpp index 0411755dda4..232ef565830 100644 --- a/components/eamxx/src/share/tests/field_tests.cpp +++ b/components/eamxx/src/share/tests/field_tests.cpp @@ -319,38 +319,31 @@ TEST_CASE("field", "") { f1.allocate_view(); randomize(f1,engine,pdf); - const int idim = 0; - const int ivar = 2; - const int sl_beg = 2; - const int sl_end = 4; - - auto f3 = f1.subfield(idim, sl_beg, sl_end); - - // const auto fl3 = f3.get_header().get_identifier().get_layout(); - // const auto alprop3 = f3.get_header().get_alloc_properties(); - // std::cout << "fl3.rank() = " << fl3.rank() << "\n"; - // auto f3dim = fl3.dims(); - // for (size_t i = 0; i < 4; i++) - // { - // std::cout << "f3dim(i) = " << f3dim[i] << "\n"; - // } - // auto f3ext = fl3.extents(); - // for (size_t i = 0; i < 4; i++) - // { - // std::cout << "f3ext(i) = " << f3ext(i) << "\n"; - // } - // std::cout << "alprop3.is_subfield() = " << alprop3.is_subfield() << "\n"; - - auto v3_h = f3.get_strided_view(true); + const int idim[4] = {0, 1, 2, 3}; + const int sl_beg[4] = {2, 3, 0, 9}; + const int sl_end[4] = {4, 6, 1, 15}; + auto v4d_h = f1.get_view(); - for (size_t i = sl_beg; i < sl_end; i++) { - for (size_t j = 0; j < d1[1]; j++) { - for (size_t k = 0; k < d1[2]; k++) { - for (size_t l = 0; l < d1[3]; l++) { - // std::cout << "v4d_h(i, j, k, l) = " << v4d_h(i, j, k, l) << "\n"; - // std::cout << "v3_h(i - sl_beg, j, k, l) = " << v3_h(i - sl_beg, j, k, l) << "\n"; - REQUIRE(v4d_h(i, j, k, l) == v3_h(i - sl_beg, j, k, l)); + int i1, i2, j1, j2, k1, k2, l1, l2; + + for (int ens = 0; ens < 4; ens++) { + auto sf = f1.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + auto sv_h = sf.get_multi_sliced_view(); + i1 = (ens == 0) ? sl_beg[0] : 0; + i2 = (ens == 0) ? sl_end[0] : d1[0]; + j1 = (ens == 1) ? sl_beg[1] : 0; + j2 = (ens == 1) ? sl_end[1] : d1[1]; + k1 = (ens == 2) ? sl_beg[2] : 0; + k2 = (ens == 2) ? sl_end[2] : d1[2]; + l1 = (ens == 3) ? sl_beg[3] : 0; + l2 = (ens == 3) ? sl_end[3] : d1[3]; + for (int i = i1; i < i2; i++) { + for (int j = j1; j < j2; j++) { + for (int k = k1; k < k2; k++) { + for (int l = l1; l < l2; l++) { + REQUIRE(v4d_h(i, j, k, l) == sv_h(i - i1, j - j1, k - k1, l - l1)); + } } } } From bcda3c0dc2ac444eb3a565d261941f392827a1a3 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Mon, 8 Apr 2024 18:23:41 -0600 Subject: [PATCH 04/16] tests finished for multi-slice sub{field,view} --- .../eamxx/src/share/field/field_impl.hpp | 44 +- .../eamxx/src/share/tests/CMakeLists.txt | 3 + .../eamxx/src/share/tests/field_tests.cpp | 149 ------- .../eamxx/src/share/tests/subfield_tests.cpp | 401 ++++++++++++++++++ 4 files changed, 426 insertions(+), 171 deletions(-) create mode 100644 components/eamxx/src/share/tests/subfield_tests.cpp diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index b938cb6070f..b0854569b24 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -204,6 +204,8 @@ auto Field::get_strided_view () const } // NOTE: multi-slicing a view is only supported for strided view return type +// TODO: N doesn't actually need to be passed in, since the subview must have +// the same rank as its parent template auto Field::get_multi_sliced_view () const -> get_strided_view_type, HD> @@ -218,9 +220,6 @@ auto Field::get_multi_sliced_view () const const auto& fl = m_header->get_identifier().get_layout(); // Checks - // TODO: decide whether a dynamic-rank view is ok - // EKAT_REQUIRE_MSG (DstRankDynamic == 1, - // "Error! Strided view not allowed with compile-time dimensions.\n"); EKAT_REQUIRE_MSG (N == fl.rank(), "Error! Input Rank must be equal to parent view's rank for multi-sliced subview.\n"); EKAT_REQUIRE_MSG(is_allocated(), @@ -230,26 +229,27 @@ auto Field::get_multi_sliced_view () const EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), "Error! Source field allocation is not compatible with the requested value type.\n"); - // Check if this field is a subview of another field + // get parent header and check if this field is a subview of another field const auto parent = m_header->get_parent().lock(); - if (parent != nullptr) { - Field f; - f.m_header = parent; - f.m_data = m_data; - - auto v_fullsize = f.get_ND_view(); - - // Now we can subview v_np1 at the correct slice - const auto& info = m_header->get_alloc_properties().get_subview_info(); - const int idim = info.dim_idx; - const int k = info.slice_idx; - const int k_end = info.slice_idx_end; - - // this version of ekat::subview overloaded conveniently, so only the single - // version of the call is required here - return DstView(ekat::subview(v_fullsize, - Kokkos::make_pair(k, k_end), idim)); - } + EKAT_REQUIRE_MSG(parent != nullptr, + "Error! Multi-sliced subview is unavailable for non-subfields.\n"); + Field f; + // create new field with the header/data from the parent + f.m_header = parent; + f.m_data = m_data; + + auto v_fullsize = f.get_ND_view(); + + // Now we can subview v_np1 at the correct slice + const auto& info = m_header->get_alloc_properties().get_subview_info(); + const int idim = info.dim_idx; + const int k_beg = info.slice_idx; + const int k_end = info.slice_idx_end; + + // this version of ekat::subview overloaded conveniently, so only the single + // version of the call is required here + return DstView(ekat::subview(v_fullsize, + Kokkos::make_pair(k_beg, k_end), idim)); } template diff --git a/components/eamxx/src/share/tests/CMakeLists.txt b/components/eamxx/src/share/tests/CMakeLists.txt index bfdc0bba2b7..4aab3ef95da 100644 --- a/components/eamxx/src/share/tests/CMakeLists.txt +++ b/components/eamxx/src/share/tests/CMakeLists.txt @@ -14,6 +14,9 @@ if (NOT SCREAM_ONLY_GENERATE_BASELINES) # Test fields CreateUnitTest(field "field_tests.cpp") + # Test fields + CreateUnitTest(subfield "subfield_tests.cpp") + # Test field utils CreateUnitTest(field_utils "field_utils.cpp" MPI_RANKS 1 ${SCREAM_TEST_MAX_RANKS}) diff --git a/components/eamxx/src/share/tests/field_tests.cpp b/components/eamxx/src/share/tests/field_tests.cpp index 232ef565830..1cb69f474de 100644 --- a/components/eamxx/src/share/tests/field_tests.cpp +++ b/components/eamxx/src/share/tests/field_tests.cpp @@ -276,155 +276,6 @@ TEST_CASE("field", "") { } } - // Subfields - SECTION ("subfield") { - std::vector t1 = {COL,CMP,CMP,LEV}; - std::vector d1 = {3,10,2,24}; - - FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); - - Field f1(fid1); - f1.allocate_view(); - randomize(f1,engine,pdf); - - const int idim = 1; - const int ivar = 2; - - auto f2 = f1.subfield(idim,ivar); - - // Wrong rank for the subfield f2 - REQUIRE_THROWS(f2.get_view()); - - auto v4d_h = f1.get_view(); - auto v3d_h = f2.get_view(); - for (int i=0; i t1 = {COL,CMP,CMP,LEV}; - std::vector d1 = {5,10,2,24}; - - FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); - - Field f1(fid1); - f1.allocate_view(); - randomize(f1,engine,pdf); - - const int idim[4] = {0, 1, 2, 3}; - const int sl_beg[4] = {2, 3, 0, 9}; - const int sl_end[4] = {4, 6, 1, 15}; - - auto v4d_h = f1.get_view(); - - int i1, i2, j1, j2, k1, k2, l1, l2; - - for (int ens = 0; ens < 4; ens++) { - auto sf = f1.subfield(idim[ens], sl_beg[ens], sl_end[ens]); - auto sv_h = sf.get_multi_sliced_view(); - i1 = (ens == 0) ? sl_beg[0] : 0; - i2 = (ens == 0) ? sl_end[0] : d1[0]; - j1 = (ens == 1) ? sl_beg[1] : 0; - j2 = (ens == 1) ? sl_end[1] : d1[1]; - k1 = (ens == 2) ? sl_beg[2] : 0; - k2 = (ens == 2) ? sl_end[2] : d1[2]; - l1 = (ens == 3) ? sl_beg[3] : 0; - l2 = (ens == 3) ? sl_end[3] : d1[3]; - for (int i = i1; i < i2; i++) { - for (int j = j1; j < j2; j++) { - for (int k = k1; k < k2; k++) { - for (int l = l1; l < l2; l++) { - REQUIRE(v4d_h(i, j, k, l) == sv_h(i - i1, j - j1, k - k1, l - l1)); - } - } - } - } - } - } - - // Dynamic Subfields - SECTION ("dynamic_subfield") { - const int vec_dim = 10; - std::vector t1 = {COL,CMP,CMP,LEV}; - std::vector d1 = {3,vec_dim,2,24}; - - FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); - - Field f1(fid1); - f1.allocate_view(); - randomize(f1,engine,pdf); - - const int idim = 1; - const int ivar = 0; - - auto f2 = f1.subfield(idim,ivar,/* dynamic = */ true); - - // Cannot reset subview idx of non-subfield fields - REQUIRE_THROWS(f1.get_header().get_alloc_properties().reset_subview_idx(0)); - - // subview idx out of bounds - auto& f2_ap = f2.get_header().get_alloc_properties(); - REQUIRE_THROWS(f2_ap.reset_subview_idx(-1)); - REQUIRE_THROWS(f2_ap.reset_subview_idx(vec_dim)); - - // Fill f1 with random numbers, and verify corresponding subviews get same values - randomize(f1,engine,pdf); - - for (int ivar_dyn=0; ivar_dyn(); - auto v3d_h = f2.get_view(); - for (int i=0; i tags_2 = {COL,CMP,LEV}; - std::vector dims_2 = {3,2,24}; - - FieldIdentifier fid_2("vec_3d",{tags_2,dims_2},m/s,"some_grid"); - - Field f_vec(fid_2); - f_vec.allocate_view(); - - auto f0 = f_vec.get_component(0); - auto f1 = f_vec.get_component(1); - - // No 3rd component - REQUIRE_THROWS(f_vec.get_component(2)); - - // f0 is scalar, no vector dimension - REQUIRE_THROWS(f0.get_component(0)); - - f0.deep_copy(1.0); - f1.deep_copy(2.0); - - f_vec.sync_to_host(); - - auto v = f_vec.get_view(); - for (int col=0; col<3; ++col) { - for (int lev=0; lev<24; ++lev) { - REQUIRE (v(col,0,lev)==1.0); - REQUIRE (v(col,1,lev)==2.0); - } - } - } - SECTION ("host_view") { Field f(fid); diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp new file mode 100644 index 00000000000..28e7ddf6129 --- /dev/null +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -0,0 +1,401 @@ +#include +#include + +#include "ekat/kokkos/ekat_subview_utils.hpp" +#include "share/field/field_identifier.hpp" +#include "share/field/field_header.hpp" +#include "share/field/field.hpp" +#include "share/field/field_manager.hpp" +#include "share/field/field_utils.hpp" +#include "share/util/scream_setup_random_test.hpp" + +#include "share/grid/point_grid.hpp" + +#include "ekat/ekat_pack.hpp" +#include "ekat/ekat_pack_utils.hpp" +#include "ekat/util/ekat_test_utils.hpp" + +namespace { + + +TEST_CASE("field", "") { + using namespace scream; + using namespace ShortFieldTagsNames; + using namespace ekat::units; + + auto engine = setup_random_test (); + using RPDF = std::uniform_real_distribution; + RPDF pdf(0.01,0.99); + + std::vector tags = {COL,LEV}; + std::vector dims = {3,24}; + + FieldIdentifier fid ("field_1", {tags,dims}, m/s,"some_grid"); + + // Subfields + SECTION ("subfield") { + std::vector t1 = {COL,CMP,CMP,LEV}; + std::vector d1 = {3,10,2,24}; + + FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); + + Field f1(fid1); + f1.allocate_view(); + randomize(f1,engine,pdf); + + const int idim = 1; + const int ivar = 2; + + auto f2 = f1.subfield(idim,ivar); + + // Wrong rank for the subfield f2 + REQUIRE_THROWS(f2.get_view()); + + auto v4d_h = f1.get_view(); + auto v3d_h = f2.get_view(); + for (int i=0; i t1 = {COL}; + std::vector d1 = {11}; + FieldIdentifier fid1("1d", {t1, d1}, m/s, "some_grid"); + + Field f1(fid1); + f1.allocate_view(); + randomize(f1, engine, pdf); + + const int idim = {0}; + const int sl_beg = {3}; + const int sl_end = {7}; + + auto v1d_h = f1.get_view(); + + // error on idx_end > v.extent(idim) + REQUIRE_THROWS(f1.subfield(idim, sl_beg, d1[0] + 1)); + // error on idx_beg negative + REQUIRE_THROWS(f1.subfield(idim, -1, sl_end)); + // error on idim > v.rank() + REQUIRE_THROWS(f1.subfield(1, sl_beg, sl_end)); + + auto sf = f1.subfield(idim, sl_beg, sl_end); + + REQUIRE_THROWS(f1.get_view()); + + auto sv_h = sf.get_multi_sliced_view(); + for (int i = sl_beg; i < sl_end; i++) { + REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); + } + } + SECTION ("2D multi-slice") { + std::vector t2 = {COL, CMP}; + std::vector d2 = {5, 10}; + FieldIdentifier fid2("2d", {t2, d2}, m/s, "some_grid"); + + Field f2(fid2); + f2.allocate_view(); + randomize(f2, engine, pdf); + + const int idim[2] = {0, 1}; + const int sl_beg[2] = {0, 3}; + const int sl_end[2] = {4, 6}; + + auto v2d_h = f2.get_view(); + int i1, i2, j1, j2; + + for (int ens = 0; ens < 2; ens++) { + auto sf = f2.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + auto sv_h = sf.get_multi_sliced_view(); + i1 = (ens == 0) ? sl_beg[0] : 0; + i2 = (ens == 0) ? sl_end[0] : d2[0]; + j1 = (ens == 1) ? sl_beg[1] : 0; + j2 = (ens == 1) ? sl_end[1] : d2[1]; + for (int i = i1; i < i2; i++) { + for (int j = j1; j < j2; j++) { + REQUIRE(v2d_h(i, j) == sv_h(i - i1, j - j1)); + } + } + } + } + SECTION ("3D multi-slice") { + // ============== + /* Rank-3 view */ + // ============== + std::vector t3 = {COL, CMP, LEV}; + std::vector d3 = {5, 10, 2}; + FieldIdentifier fid3("3d", {t3, d3}, m/s, "some_grid"); + + Field f3(fid3); + f3.allocate_view(); + randomize(f3, engine, pdf); + + const int idim[3] = {0, 1, 2}; + const int sl_beg[3] = {2, 3, 0}; + const int sl_end[3] = {4, 6, 1}; + + auto v3d_h = f3.get_view(); + int i1, i2, j1, j2, k1, k2; + + for (int ens = 0; ens < 3; ens++) { + auto sf = f3.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + auto sv_h = sf.get_multi_sliced_view(); + i1 = (ens == 0) ? sl_beg[0] : 0; + i2 = (ens == 0) ? sl_end[0] : d3[0]; + j1 = (ens == 1) ? sl_beg[1] : 0; + j2 = (ens == 1) ? sl_end[1] : d3[1]; + k1 = (ens == 2) ? sl_beg[2] : 0; + k2 = (ens == 2) ? sl_end[2] : d3[2]; + for (int i = i1; i < i2; i++) { + for (int j = j1; j < j2; j++) { + for (int k = k1; k < k2; k++) { + REQUIRE(v3d_h(i, j, k) == sv_h(i - i1, j - j1, k - k1)); + } + } + } + } + // ====================== + /* get_component test */ + // ====================== + auto cmp3 = f3.get_component(sl_beg[1], sl_end[1]); + REQUIRE_THROWS(f3.get_component(sl_beg[1], d3[1] + 1)); + REQUIRE_THROWS(f3.get_component(-1, sl_end[1] + 1)); + + auto svc_h = cmp3.get_multi_sliced_view(); + for (int i = 0; i < d3[0]; i++) { + for (int j = sl_beg[1]; j < sl_end[1]; j++) { + for (int k = 0; k < d3[2]; k++) { + REQUIRE(v3d_h(i, j, k) == svc_h(i, j - sl_beg[1], k)); + } + } + } + + } + SECTION ("4D multi-slice") { + // ============== + /* Rank-4 view */ + // ============== + std::vector t4 = {COL, CMP, CMP, LEV}; + std::vector d4 = {5, 10, 2, 23}; + FieldIdentifier fid4("4d", {t4, d4}, m/s, "some_grid"); + + Field f4(fid4); + f4.allocate_view(); + randomize(f4, engine, pdf); + + const int idim[4] = {0, 1, 2, 3}; + const int sl_beg[4] = {2, 3, 0, 9}; + const int sl_end[4] = {4, 6, 1, 15}; + + auto v4d_h = f4.get_view(); + int i1, i2, j1, j2, k1, k2, l1, l2; + + for (int ens = 0; ens < 4; ens++) { + auto sf = f4.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + auto sv_h = sf.get_multi_sliced_view(); + i1 = (ens == 0) ? sl_beg[0] : 0; + i2 = (ens == 0) ? sl_end[0] : d4[0]; + j1 = (ens == 1) ? sl_beg[1] : 0; + j2 = (ens == 1) ? sl_end[1] : d4[1]; + k1 = (ens == 2) ? sl_beg[2] : 0; + k2 = (ens == 2) ? sl_end[2] : d4[2]; + l1 = (ens == 3) ? sl_beg[3] : 0; + l2 = (ens == 3) ? sl_end[3] : d4[3]; + for (int i = i1; i < i2; i++) { + for (int j = j1; j < j2; j++) { + for (int k = k1; k < k2; k++) { + for (int l = l1; l < l2; l++) { + REQUIRE(v4d_h(i, j, k, l) == sv_h(i - i1, j - j1, k - k1, l - l1)); + } + } + } + } + } + } + SECTION ("5D multi-slice") { + // ============== + /* Rank-5 view */ + // ============== + std::vector t5 = {EL, CMP, GP, GP, LEV}; + std::vector d5 = {5, 10, 4, 2, 23}; + FieldIdentifier fid5("5d", {t5, d5}, m/s, "some_grid"); + + Field f5(fid5); + f5.allocate_view(); + randomize(f5, engine, pdf); + + const int idim[5] = {0, 1, 2, 3, 4}; + const int sl_beg[5] = {2, 3, 1, 0, 9}; + const int sl_end[5] = {4, 6, 3, 1, 15}; + + auto v5d_h = f5.get_view(); + int i1, i2, j1, j2, k1, k2, l1, l2, m1, m2; + + for (int ens = 0; ens < 5; ens++) { + auto sf = f5.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + auto sv_h = sf.get_multi_sliced_view(); + i1 = (ens == 0) ? sl_beg[0] : 0; + i2 = (ens == 0) ? sl_end[0] : d5[0]; + j1 = (ens == 1) ? sl_beg[1] : 0; + j2 = (ens == 1) ? sl_end[1] : d5[1]; + k1 = (ens == 2) ? sl_beg[2] : 0; + k2 = (ens == 2) ? sl_end[2] : d5[2]; + l1 = (ens == 3) ? sl_beg[3] : 0; + l2 = (ens == 3) ? sl_end[3] : d5[3]; + m1 = (ens == 4) ? sl_beg[4] : 0; + m2 = (ens == 4) ? sl_end[4] : d5[4]; + for (int i = i1; i < i2; i++) { + for (int j = j1; j < j2; j++) { + for (int k = k1; k < k2; k++) { + for (int l = l1; l < l2; l++) { + for (int m = m1; m < m2; m++) { + REQUIRE(v5d_h(i, j, k, l, m) == sv_h(i - i1, j - j1, k - k1, + l - l1, m - m1)); + } + } + } + } + } + } + } + SECTION ("6D multi-slice") { + // ============== + /* Rank-6 view */ + // ============== + std::vector t6 = {EL, TL, CMP, GP, GP, LEV}; + std::vector d6 = {5, 10, 4, 2, 9, 23}; + FieldIdentifier fid6("6d", {t6, d6}, m/s, "some_grid"); + + Field f6(fid6); + f6.allocate_view(); + randomize(f6, engine, pdf); + + const int idim[6] = {0, 1, 2, 3, 4, 5}; + const int sl_beg[6] = {2, 3, 1, 0, 5, 9}; + const int sl_end[6] = {4, 6, 3, 1, 8, 15}; + + auto v6d_h = f6.get_view(); + int i1, i2, j1, j2, k1, k2, l1, l2, m1, m2, n1, n2; + + for (int ens = 0; ens < 6; ens++) { + auto sf = f6.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + auto sv_h = sf.get_multi_sliced_view(); + i1 = (ens == 0) ? sl_beg[0] : 0; + i2 = (ens == 0) ? sl_end[0] : d6[0]; + j1 = (ens == 1) ? sl_beg[1] : 0; + j2 = (ens == 1) ? sl_end[1] : d6[1]; + k1 = (ens == 2) ? sl_beg[2] : 0; + k2 = (ens == 2) ? sl_end[2] : d6[2]; + l1 = (ens == 3) ? sl_beg[3] : 0; + l2 = (ens == 3) ? sl_end[3] : d6[3]; + m1 = (ens == 4) ? sl_beg[4] : 0; + m2 = (ens == 4) ? sl_end[4] : d6[4]; + n1 = (ens == 5) ? sl_beg[5] : 0; + n2 = (ens == 5) ? sl_end[5] : d6[5]; + for (int i = i1; i < i2; i++) { + for (int j = j1; j < j2; j++) { + for (int k = k1; k < k2; k++) { + for (int l = l1; l < l2; l++) { + for (int m = m1; m < m2; m++) { + for (int n = n1; n < n2; n++) { + REQUIRE(v6d_h(i, j, k, l, m, n) == + sv_h(i - i1, j - j1, k - k1, + l - l1, m - m1, n - n1)); + } + } + } + } + } + } + } + } + } + + // Dynamic Subfields + SECTION ("dynamic_subfield") { + const int vec_dim = 10; + std::vector t1 = {COL,CMP,CMP,LEV}; + std::vector d1 = {3,vec_dim,2,24}; + + FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); + + Field f1(fid1); + f1.allocate_view(); + randomize(f1,engine,pdf); + + const int idim = 1; + const int ivar = 0; + + auto f2 = f1.subfield(idim,ivar,/* dynamic = */ true); + + // Cannot reset subview idx of non-subfield fields + REQUIRE_THROWS(f1.get_header().get_alloc_properties().reset_subview_idx(0)); + + // subview idx out of bounds + auto& f2_ap = f2.get_header().get_alloc_properties(); + REQUIRE_THROWS(f2_ap.reset_subview_idx(-1)); + REQUIRE_THROWS(f2_ap.reset_subview_idx(vec_dim)); + + // Fill f1 with random numbers, and verify corresponding subviews get same values + randomize(f1,engine,pdf); + + for (int ivar_dyn=0; ivar_dyn(); + auto v3d_h = f2.get_view(); + for (int i=0; i tags_2 = {COL,CMP,LEV}; + std::vector dims_2 = {3,2,24}; + + FieldIdentifier fid_2("vec_3d",{tags_2,dims_2},m/s,"some_grid"); + + Field f_vec(fid_2); + f_vec.allocate_view(); + + auto f0 = f_vec.get_component(0); + auto f1 = f_vec.get_component(1); + + // No 3rd component + REQUIRE_THROWS(f_vec.get_component(2)); + + // f0 is scalar, no vector dimension + REQUIRE_THROWS(f0.get_component(0)); + + f0.deep_copy(1.0); + f1.deep_copy(2.0); + + f_vec.sync_to_host(); + + auto v = f_vec.get_view(); + for (int col=0; col<3; ++col) { + for (int lev=0; lev<24; ++lev) { + REQUIRE (v(col,0,lev)==1.0); + REQUIRE (v(col,1,lev)==2.0); + } + } + } +} + + + + + +} // anonymous namespace From b1a19c1bdbd9d117b19f09e272507f20e39bb996 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Wed, 10 Apr 2024 22:16:06 -0600 Subject: [PATCH 05/16] address PR comments --- .../eamxx/src/share/field/field_impl.hpp | 39 ++++++++++--------- .../eamxx/src/share/tests/subfield_tests.cpp | 2 +- 2 files changed, 21 insertions(+), 20 deletions(-) diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index b0854569b24..91c196bc465 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -206,10 +206,9 @@ auto Field::get_strided_view () const // NOTE: multi-slicing a view is only supported for strided view return type // TODO: N doesn't actually need to be passed in, since the subview must have // the same rank as its parent -template -auto Field::get_multi_sliced_view () const - -> get_strided_view_type, HD> -{ +template +auto Field::get_multi_sliced_view() const + -> get_strided_view_type, HD> { // The destination view type on correct mem space using DstView = get_strided_view_type, HD>; // The dst value types @@ -220,18 +219,22 @@ auto Field::get_multi_sliced_view () const const auto& fl = m_header->get_identifier().get_layout(); // Checks - EKAT_REQUIRE_MSG (N == fl.rank(), - "Error! Input Rank must be equal to parent view's rank for multi-sliced subview.\n"); - EKAT_REQUIRE_MSG(is_allocated(), + EKAT_REQUIRE_MSG(N == fl.rank(), "Error! Input Rank must be equal to parent " + "view's rank for multi-sliced subview.\n"); + EKAT_REQUIRE_MSG( + is_allocated(), "Error! Cannot extract a field's view before allocation happens.\n"); - EKAT_REQUIRE_MSG (not m_is_read_only || std::is_const::value, - "Error! Cannot get a view to non-const data if the field is read-only.\n"); + EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, + "Error! Cannot get a view to non-const data if the field is " + "read-only.\n"); EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), - "Error! Source field allocation is not compatible with the requested value type.\n"); + "Error! Source field allocation is not compatible with the " + "requested value type.\n"); // get parent header and check if this field is a subview of another field const auto parent = m_header->get_parent().lock(); - EKAT_REQUIRE_MSG(parent != nullptr, + EKAT_REQUIRE_MSG( + parent != nullptr, "Error! Multi-sliced subview is unavailable for non-subfields.\n"); Field f; // create new field with the header/data from the parent @@ -242,14 +245,12 @@ auto Field::get_multi_sliced_view () const // Now we can subview v_np1 at the correct slice const auto& info = m_header->get_alloc_properties().get_subview_info(); - const int idim = info.dim_idx; - const int k_beg = info.slice_idx; - const int k_end = info.slice_idx_end; - - // this version of ekat::subview overloaded conveniently, so only the single - // version of the call is required here - return DstView(ekat::subview(v_fullsize, - Kokkos::make_pair(k_beg, k_end), idim)); + const int idim = info.dim_idx; + const int k_beg = info.slice_idx; + const int k_end = info.slice_idx_end; + + return DstView(ekat::subview( + v_fullsize, Kokkos::make_pair(k_beg, k_end), idim)); } template diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp index 28e7ddf6129..24950388fa5 100644 --- a/components/eamxx/src/share/tests/subfield_tests.cpp +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -88,7 +88,7 @@ TEST_CASE("field", "") { auto sf = f1.subfield(idim, sl_beg, sl_end); - REQUIRE_THROWS(f1.get_view()); + REQUIRE_THROWS(f1.get_multi_sliced_view()); auto sv_h = sf.get_multi_sliced_view(); for (int i = sl_beg; i < sl_end; i++) { From e57794738560417c4f8566dfba960e5a91a63e6b Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Thu, 11 Apr 2024 15:57:02 -0600 Subject: [PATCH 06/16] remove require_throws() from subfield tests, and point ekat to master after merge --- .../eamxx/src/share/tests/subfield_tests.cpp | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp index 24950388fa5..fd4203feefb 100644 --- a/components/eamxx/src/share/tests/subfield_tests.cpp +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -78,19 +78,9 @@ TEST_CASE("field", "") { const int sl_end = {7}; auto v1d_h = f1.get_view(); - - // error on idx_end > v.extent(idim) - REQUIRE_THROWS(f1.subfield(idim, sl_beg, d1[0] + 1)); - // error on idx_beg negative - REQUIRE_THROWS(f1.subfield(idim, -1, sl_end)); - // error on idim > v.rank() - REQUIRE_THROWS(f1.subfield(1, sl_beg, sl_end)); - - auto sf = f1.subfield(idim, sl_beg, sl_end); - - REQUIRE_THROWS(f1.get_multi_sliced_view()); - + auto sf = f1.subfield(idim, sl_beg, sl_end auto sv_h = sf.get_multi_sliced_view(); + for (int i = sl_beg; i < sl_end; i++) { REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); } @@ -165,8 +155,6 @@ TEST_CASE("field", "") { /* get_component test */ // ====================== auto cmp3 = f3.get_component(sl_beg[1], sl_end[1]); - REQUIRE_THROWS(f3.get_component(sl_beg[1], d3[1] + 1)); - REQUIRE_THROWS(f3.get_component(-1, sl_end[1] + 1)); auto svc_h = cmp3.get_multi_sliced_view(); for (int i = 0; i < d3[0]; i++) { From f7ee547c3048e37bde63b1cdc351f1bac0d7d3e1 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Tue, 16 Apr 2024 15:33:46 -0600 Subject: [PATCH 07/16] fix typo --- components/eamxx/src/share/tests/subfield_tests.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp index fd4203feefb..5c649f9a272 100644 --- a/components/eamxx/src/share/tests/subfield_tests.cpp +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -78,7 +78,7 @@ TEST_CASE("field", "") { const int sl_end = {7}; auto v1d_h = f1.get_view(); - auto sf = f1.subfield(idim, sl_beg, sl_end + auto sf = f1.subfield(idim, sl_beg, sl_end); auto sv_h = sf.get_multi_sliced_view(); for (int i = sl_beg; i < sl_end; i++) { From 6d14e3a90e170a5fdc1077c084731b20e9b4f6bb Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Wed, 17 Apr 2024 15:13:01 -0600 Subject: [PATCH 08/16] edits in response to Luca review --- components/eamxx/src/share/field/field.cpp | 27 ++- components/eamxx/src/share/field/field.hpp | 13 +- .../src/share/field/field_alloc_prop.cpp | 15 +- .../src/share/field/field_alloc_prop.hpp | 28 ++-- .../eamxx/src/share/field/field_header.cpp | 4 +- .../eamxx/src/share/field/field_header.hpp | 6 +- .../eamxx/src/share/tests/subfield_tests.cpp | 154 ++++++++---------- 7 files changed, 112 insertions(+), 135 deletions(-) diff --git a/components/eamxx/src/share/field/field.cpp b/components/eamxx/src/share/field/field.cpp index 1e790ae3cea..f56b3b83d93 100644 --- a/components/eamxx/src/share/field/field.cpp +++ b/components/eamxx/src/share/field/field.cpp @@ -82,7 +82,6 @@ Field Field:: subfield (const std::string& sf_name, const ekat::units::Units& sf_units, const int idim, const int index, const bool dynamic) const { - // const auto& id = m_header->get_identifier(); const auto& lt = id.get_layout(); @@ -119,8 +118,7 @@ subfield (const int idim, const int index, const bool dynamic) const { // slice at index idim, entries \in [index_beg, index_end] Field Field::subfield(const std::string& sf_name, const ekat::units::Units& sf_units, const int idim, - const int index_beg, const int index_end, - const bool dynamic) const { + const int index_beg, const int index_end) const { const auto& id = m_header->get_identifier(); const auto& lt = id.get_layout(); @@ -129,11 +127,9 @@ Field Field::subfield(const std::string& sf_name, EKAT_REQUIRE_MSG( is_allocated(), "Error! Input field must be allocated in order to subview it.\n"); - // EKAT_REQUIRE_MSG(idim == 0 || idim == 1, - // "Error! Subview dimension index must be either 0 or 1.\n"); auto sf_layout = - lt.clone_with_different_extent(idim, index_end - index_beg + 1); + lt.clone_with_different_extent(idim, index_end - index_beg); // Create identifier for subfield FieldIdentifier sf_id(sf_name, sf_layout, sf_units, id.get_grid_name()); @@ -141,23 +137,22 @@ Field Field::subfield(const std::string& sf_name, // Note: we can access protected members, since it's the same type Field sf; sf.m_header = create_subfield_header(sf_id, m_header, idim, index_beg, - index_end, dynamic); + index_end); sf.m_data = m_data; return sf; } Field Field::subfield(const std::string& sf_name, const int idim, - const int index_beg, const int index_end, - const bool dynamic) const { + const int index_beg, const int index_end) const { const auto& id = m_header->get_identifier(); - return subfield(sf_name, id.get_units(), idim, index_beg, index_end, dynamic); + return subfield(sf_name, id.get_units(), idim, index_beg, index_end); } -Field Field::subfield(const int idim, const int index_beg, const int index_end, - const bool dynamic) const { - return subfield(m_header->get_identifier().name(), idim, index_beg, index_end, - dynamic); +Field Field::subfield(const int idim, const int index_beg, + const int index_end) const { + return subfield(m_header->get_identifier().name(), idim, index_beg, + index_end); } Field Field:: @@ -177,7 +172,7 @@ get_component (const int i, const bool dynamic) { return subfield (fname + "_" + std::to_string(i),idim,i,dynamic); } -Field Field::get_component(const int i1, const int i2, const bool dynamic) { +Field Field::get_components(const int i1, const int i2) { const auto& layout = get_header().get_identifier().get_layout(); const auto& fname = get_header().get_identifier().name(); EKAT_REQUIRE_MSG(layout.is_vector_layout(), @@ -195,7 +190,7 @@ Field Field::get_component(const int i1, const int i2, const bool dynamic) { // Add _$i1-$i2 to the field name, to avoid issues if the subfield is stored // in some structure that requires unique names (e.g., a remapper) return subfield(fname + "_" + std::to_string(i1) + "-" + std::to_string(i2), - idim, i1, i2, dynamic); + idim, i1, i2); } bool Field::equivalent(const Field& rhs) const diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index 98169ece28d..f4ad8199726 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -247,23 +247,20 @@ class Field { Field subfield (const std::string& sf_name, const int idim, const int index, const bool dynamic = false) const; Field subfield (const int idim, const int k, const bool dynamic = false) const; - // subfield fxn to extract multiple slices in a continuous range of indices + // extracts a subfield composed of multiple slices in a continuous range of indices // e.g., (in matlab syntax) subf = f.subfield(:, 1:3, :) // but NOT subf = f.subfield(:, [1, 3, 4], :) Field subfield (const std::string& sf_name, const ekat::units::Units& sf_units, - const int idim, const int index_beg, const int index_end, - const bool dynamic = false) const; + const int idim, const int index_beg, const int index_end) const; Field subfield (const std::string& sf_name, const int idim, - const int index_beg, const int index_end, - const bool dynamic = false) const; - Field subfield (const int idim, const int k_beg, const int k_end, - const bool dynamic = false) const; + const int index_beg, const int index_end) const; + Field subfield (const int idim, const int k_beg, const int k_end) const; // If this field is a vector field, get a subfield for the ith component. // If dynamic = true, it is possible to "reset" the component index at runtime. // Note: throws if this is not a vector field. Field get_component (const int i, const bool dynamic = false); // version for slicing across multiple, contiguous indices - Field get_component (const int i1, const int i2, const bool dynamic = false); + Field get_components (const int i1, const int i2); // Checks whether the underlying view has been already allocated. bool is_allocated () const { return m_data.d_view.data()!=nullptr; } diff --git a/components/eamxx/src/share/field/field_alloc_prop.cpp b/components/eamxx/src/share/field/field_alloc_prop.cpp index 08d850a9325..635e62bae19 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.cpp +++ b/components/eamxx/src/share/field/field_alloc_prop.cpp @@ -79,12 +79,11 @@ subview (const int idim, const int k, const bool dynamic) const { FieldAllocProp FieldAllocProp::subview(const int idim, const int k_beg, - const int k_end, - const bool dynamic) const { + const int k_end) const { EKAT_REQUIRE_MSG( is_committed(), "Error! Subview requires alloc properties to be committed.\n"); - EKAT_REQUIRE_MSG(idim <= m_layout.rank(), + EKAT_REQUIRE_MSG(idim < m_layout.rank(), "Error! Dimension index out of bounds.\n"); EKAT_REQUIRE_MSG(k_beg < k_end, "Error! Slice indices are invalid (non-increasing).\n"); @@ -97,19 +96,17 @@ FieldAllocProp FieldAllocProp::subview(const int idim, props.m_committed = true; props.m_scalar_type_size = m_scalar_type_size; props.m_layout = - m_layout.clone_with_different_extent(idim, k_end - k_beg + 1); + m_layout.clone_with_different_extent(idim, k_end - k_beg); // Output is contiguous if either // - this->m_contiguous=true AND idim==0 // - m_layout.dim(i)==1 for all i create_subfield_header(const FieldIdentifier& id, std::shared_ptr parent, const int idim, - const int k_beg, const int k_end, const bool dynamic) { + const int k_beg, const int k_end) { // Sanity checks EKAT_REQUIRE_MSG(parent != nullptr, "Error! Invalid pointer for parent header.\n"); @@ -92,7 +92,7 @@ create_subfield_header(const FieldIdentifier& id, // Create alloc props fh->m_alloc_prop = std::make_shared( - parent->get_alloc_properties().subview(idim, k_beg, k_end, dynamic)); + parent->get_alloc_properties().subview(idim, k_beg, k_end)); return fh; } diff --git a/components/eamxx/src/share/field/field_header.hpp b/components/eamxx/src/share/field/field_header.hpp index 84284272013..fce98965e2f 100644 --- a/components/eamxx/src/share/field/field_header.hpp +++ b/components/eamxx/src/share/field/field_header.hpp @@ -95,8 +95,7 @@ class FieldHeader : public FamilyTracking { friend std::shared_ptr create_subfield_header (const FieldIdentifier&, std::shared_ptr, - const int idim, const int k_beg, const int k_end, - const bool dynamic); + const int idim, const int k_beg, const int k_end); // NOTE: the identifier *cannot* be a shared_ptr, b/c we // don't foresee sharing an identifier between two @@ -178,8 +177,7 @@ create_subfield_header (const FieldIdentifier& id, std::shared_ptr create_subfield_header (const FieldIdentifier& id, std::shared_ptr parent, - const int idim, const int k_beg, const int k_end, - const bool dynamic); + const int idim, const int k_beg, const int k_end); } // namespace scream diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp index 5c649f9a272..b68e62b4fdd 100644 --- a/components/eamxx/src/share/tests/subfield_tests.cpp +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -1,73 +1,59 @@ #include #include -#include "ekat/kokkos/ekat_subview_utils.hpp" -#include "share/field/field_identifier.hpp" -#include "share/field/field_header.hpp" #include "share/field/field.hpp" -#include "share/field/field_manager.hpp" #include "share/field/field_utils.hpp" #include "share/util/scream_setup_random_test.hpp" -#include "share/grid/point_grid.hpp" - -#include "ekat/ekat_pack.hpp" -#include "ekat/ekat_pack_utils.hpp" #include "ekat/util/ekat_test_utils.hpp" namespace { - TEST_CASE("field", "") { using namespace scream; using namespace ShortFieldTagsNames; using namespace ekat::units; - auto engine = setup_random_test (); + auto engine = setup_random_test(); using RPDF = std::uniform_real_distribution; - RPDF pdf(0.01,0.99); - - std::vector tags = {COL,LEV}; - std::vector dims = {3,24}; - - FieldIdentifier fid ("field_1", {tags,dims}, m/s,"some_grid"); + RPDF pdf(0.01, 0.99); // Subfields - SECTION ("subfield") { - std::vector t1 = {COL,CMP,CMP,LEV}; - std::vector d1 = {3,10,2,24}; + SECTION("subfield") { + std::vector t1 = {COL, CMP, CMP, LEV}; + std::vector d1 = {3, 10, 2, 24}; - FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); + FieldIdentifier fid1("4d", {t1, d1}, m / s, "some_grid"); Field f1(fid1); f1.allocate_view(); - randomize(f1,engine,pdf); + randomize(f1, engine, pdf); const int idim = 1; const int ivar = 2; - auto f2 = f1.subfield(idim,ivar); + auto f2 = f1.subfield(idim, ivar); // Wrong rank for the subfield f2 REQUIRE_THROWS(f2.get_view()); - auto v4d_h = f1.get_view(); - auto v3d_h = f2.get_view(); - for (int i=0; i(); + auto v3d_h = f2.get_view(); + for (int i = 0; i < d1[0]; ++i) + for (int j = 0; j < d1[2]; ++j) + for (int k = 0; k < d1[3]; ++k) { + REQUIRE(v4d_h(i, ivar, j, k) == v3d_h(i, j, k)); } } - SECTION ("multi-sliced subfield") { - SECTION ("1D multi-slice") { + SECTION("multi-sliced subfield") { + SECTION("1D multi-slice") { // ============== /* Rank-1 view */ // ============== std::vector t1 = {COL}; std::vector d1 = {11}; - FieldIdentifier fid1("1d", {t1, d1}, m/s, "some_grid"); + FieldIdentifier fid1("1d", {t1, d1}, m / s, "some_grid"); Field f1(fid1); f1.allocate_view(); @@ -81,14 +67,17 @@ TEST_CASE("field", "") { auto sf = f1.subfield(idim, sl_beg, sl_end); auto sv_h = sf.get_multi_sliced_view(); + REQUIRE(sv_h.extent_int(idim) == (sl_end - sl_beg)); + for (int i = sl_beg; i < sl_end; i++) { REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); } } - SECTION ("2D multi-slice") { + + SECTION("2D multi-slice") { std::vector t2 = {COL, CMP}; std::vector d2 = {5, 10}; - FieldIdentifier fid2("2d", {t2, d2}, m/s, "some_grid"); + FieldIdentifier fid2("2d", {t2, d2}, m / s, "some_grid"); Field f2(fid2); f2.allocate_view(); @@ -115,13 +104,13 @@ TEST_CASE("field", "") { } } } - SECTION ("3D multi-slice") { + SECTION("3D multi-slice") { // ============== /* Rank-3 view */ // ============== std::vector t3 = {COL, CMP, LEV}; std::vector d3 = {5, 10, 2}; - FieldIdentifier fid3("3d", {t3, d3}, m/s, "some_grid"); + FieldIdentifier fid3("3d", {t3, d3}, m / s, "some_grid"); Field f3(fid3); f3.allocate_view(); @@ -152,9 +141,9 @@ TEST_CASE("field", "") { } } // ====================== - /* get_component test */ + /* get_components test */ // ====================== - auto cmp3 = f3.get_component(sl_beg[1], sl_end[1]); + auto cmp3 = f3.get_components(sl_beg[1], sl_end[1]); auto svc_h = cmp3.get_multi_sliced_view(); for (int i = 0; i < d3[0]; i++) { @@ -164,15 +153,14 @@ TEST_CASE("field", "") { } } } - } - SECTION ("4D multi-slice") { + SECTION("4D multi-slice") { // ============== /* Rank-4 view */ // ============== std::vector t4 = {COL, CMP, CMP, LEV}; std::vector d4 = {5, 10, 2, 23}; - FieldIdentifier fid4("4d", {t4, d4}, m/s, "some_grid"); + FieldIdentifier fid4("4d", {t4, d4}, m / s, "some_grid"); Field f4(fid4); f4.allocate_view(); @@ -200,20 +188,21 @@ TEST_CASE("field", "") { for (int j = j1; j < j2; j++) { for (int k = k1; k < k2; k++) { for (int l = l1; l < l2; l++) { - REQUIRE(v4d_h(i, j, k, l) == sv_h(i - i1, j - j1, k - k1, l - l1)); + REQUIRE(v4d_h(i, j, k, l) == + sv_h(i - i1, j - j1, k - k1, l - l1)); } } } } } } - SECTION ("5D multi-slice") { + SECTION("5D multi-slice") { // ============== /* Rank-5 view */ // ============== std::vector t5 = {EL, CMP, GP, GP, LEV}; std::vector d5 = {5, 10, 4, 2, 23}; - FieldIdentifier fid5("5d", {t5, d5}, m/s, "some_grid"); + FieldIdentifier fid5("5d", {t5, d5}, m / s, "some_grid"); Field f5(fid5); f5.allocate_view(); @@ -244,8 +233,8 @@ TEST_CASE("field", "") { for (int k = k1; k < k2; k++) { for (int l = l1; l < l2; l++) { for (int m = m1; m < m2; m++) { - REQUIRE(v5d_h(i, j, k, l, m) == sv_h(i - i1, j - j1, k - k1, - l - l1, m - m1)); + REQUIRE(v5d_h(i, j, k, l, m) == + sv_h(i - i1, j - j1, k - k1, l - l1, m - m1)); } } } @@ -253,13 +242,13 @@ TEST_CASE("field", "") { } } } - SECTION ("6D multi-slice") { + SECTION("6D multi-slice") { // ============== /* Rank-6 view */ // ============== std::vector t6 = {EL, TL, CMP, GP, GP, LEV}; std::vector d6 = {5, 10, 4, 2, 9, 23}; - FieldIdentifier fid6("6d", {t6, d6}, m/s, "some_grid"); + FieldIdentifier fid6("6d", {t6, d6}, m / s, "some_grid"); Field f6(fid6); f6.allocate_view(); @@ -293,9 +282,15 @@ TEST_CASE("field", "") { for (int l = l1; l < l2; l++) { for (int m = m1; m < m2; m++) { for (int n = n1; n < n2; n++) { - REQUIRE(v6d_h(i, j, k, l, m, n) == - sv_h(i - i1, j - j1, k - k1, - l - l1, m - m1, n - n1)); + REQUIRE(v6d_h(i, j, k, l, m, n) == sv_h(i - i1, j - j1, + k - k1, l - l1, + m - m1, n - n1)); + REQUIRE((sv_h.extent_int(0) == (i2 - i1) && + sv_h.extent_int(1) == (j2 - j1) && + sv_h.extent_int(2) == (k2 - k1) && + sv_h.extent_int(3) == (l2 - l1) && + sv_h.extent_int(4) == (m2 - m1) && + sv_h.extent_int(5) == (n2 - n1))); } } } @@ -305,23 +300,22 @@ TEST_CASE("field", "") { } } } - // Dynamic Subfields - SECTION ("dynamic_subfield") { + SECTION("dynamic_subfield") { const int vec_dim = 10; - std::vector t1 = {COL,CMP,CMP,LEV}; - std::vector d1 = {3,vec_dim,2,24}; + std::vector t1 = {COL, CMP, CMP, LEV}; + std::vector d1 = {3, vec_dim, 2, 24}; - FieldIdentifier fid1("4d",{t1,d1},m/s,"some_grid"); + FieldIdentifier fid1("4d", {t1, d1}, m / s, "some_grid"); Field f1(fid1); f1.allocate_view(); - randomize(f1,engine,pdf); + randomize(f1, engine, pdf); const int idim = 1; const int ivar = 0; - auto f2 = f1.subfield(idim,ivar,/* dynamic = */ true); + auto f2 = f1.subfield(idim, ivar, /* dynamic = */ true); // Cannot reset subview idx of non-subfield fields REQUIRE_THROWS(f1.get_header().get_alloc_properties().reset_subview_idx(0)); @@ -331,29 +325,30 @@ TEST_CASE("field", "") { REQUIRE_THROWS(f2_ap.reset_subview_idx(-1)); REQUIRE_THROWS(f2_ap.reset_subview_idx(vec_dim)); - // Fill f1 with random numbers, and verify corresponding subviews get same values - randomize(f1,engine,pdf); + // Fill f1 with random numbers, and verify corresponding subviews get same + // values + randomize(f1, engine, pdf); - for (int ivar_dyn=0; ivar_dyn(); - auto v3d_h = f2.get_view(); - for (int i=0; i(); + auto v3d_h = f2.get_view(); + for (int i = 0; i < d1[0]; ++i) + for (int j = 0; j < d1[2]; ++j) + for (int k = 0; k < d1[3]; ++k) { + REQUIRE(v4d_h(i, ivar_dyn, j, k) == v3d_h(i, j, k)); } } } - SECTION ("vector_component") { - std::vector tags_2 = {COL,CMP,LEV}; - std::vector dims_2 = {3,2,24}; + SECTION("vector_component") { + std::vector tags_2 = {COL, CMP, LEV}; + std::vector dims_2 = {3, 2, 24}; - FieldIdentifier fid_2("vec_3d",{tags_2,dims_2},m/s,"some_grid"); + FieldIdentifier fid_2("vec_3d", {tags_2, dims_2}, m / s, "some_grid"); Field f_vec(fid_2); f_vec.allocate_view(); @@ -372,18 +367,13 @@ TEST_CASE("field", "") { f_vec.sync_to_host(); - auto v = f_vec.get_view(); - for (int col=0; col<3; ++col) { - for (int lev=0; lev<24; ++lev) { - REQUIRE (v(col,0,lev)==1.0); - REQUIRE (v(col,1,lev)==2.0); + auto v = f_vec.get_view(); + for (int col = 0; col < 3; ++col) { + for (int lev = 0; lev < 24; ++lev) { + REQUIRE(v(col, 0, lev) == 1.0); + REQUIRE(v(col, 1, lev) == 2.0); } } } } - - - - - } // anonymous namespace From 8c9384d74085684ee8929700a1fe7b6fd0b3130b Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Wed, 24 Apr 2024 13:10:13 -0600 Subject: [PATCH 09/16] wip on unifying subview --- .../eamxx/src/share/field/field_impl.hpp | 194 ++++++++++-------- .../eamxx/src/share/tests/subfield_tests.cpp | 54 ++--- 2 files changed, 136 insertions(+), 112 deletions(-) diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index 91c196bc465..c2fdce5d5cf 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -141,117 +141,141 @@ auto Field::get_view () const return DstView(view_ND); } -template -auto Field::get_strided_view () const - -> get_strided_view_type -{ +template +auto Field::get_strided_view() const -> +get_strided_view_type { // The destination view type on correct mem space - using DstView = get_strided_view_type; + using DstView = get_strided_view_type; // The dst value types using DstValueType = typename DstView::traits::value_type; // We only allow to reshape to a view of the correct rank constexpr int DstRank = DstView::rank; - constexpr int DstRankDynamic= DstView::rank_dynamic; + constexpr int DstRankDynamic = DstView::rank_dynamic; // Get src details const auto& alloc_prop = m_header->get_alloc_properties(); const auto& fl = m_header->get_identifier().get_layout(); // Checks - EKAT_REQUIRE_MSG (DstRank==1 && fl.rank()==1, - "Error! Strided view only available for rank-1 fields.\n"); - EKAT_REQUIRE_MSG (DstRankDynamic==1, - "Error! Strided view not allowed with compile-time dimensions.\n"); - EKAT_REQUIRE_MSG(is_allocated(), + EKAT_REQUIRE_MSG( + is_allocated(), "Error! Cannot extract a field's view before allocation happens.\n"); - EKAT_REQUIRE_MSG (not m_is_read_only || std::is_const::value, - "Error! Cannot get a view to non-const data if the field is read-only.\n"); + EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, + "Error! Cannot get a view to non-const data if the field is " + "read-only.\n"); EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), - "Error! Source field allocation is not compatible with the requested value type.\n"); + "Error! Source field allocation is not compatible with the " + "requested value type.\n"); // Check if this field is a subview of another field const auto parent = m_header->get_parent().lock(); - if (parent!=nullptr) { + if (parent != nullptr) { // Parent field has correct layout to reinterpret the view into N+1-dim view - // So create the parent field on the fly, use it to get the N+1-dim view, then subview it. - // NOTE: we can set protected members, since f is the same type of this class. + // So create the parent field on the fly, use it to get the N+1-dim view, + // then subview it. NOTE: we can set protected members, since f is the same + // type of this class. Field f; f.m_header = parent; - f.m_data = m_data; - - // Take 2 dimensional view with normal LayoutRight - auto v_np1 = f.get_ND_view(); + f.m_data = m_data; - // Now we can subview v_np1 at the correct slice - const auto& info = m_header->get_alloc_properties().get_subview_info(); + // get subview info to determine whether we are single- or multi-slicing + const auto& sv_alloc_prop = m_header->get_alloc_properties(); + const auto& info = sv_alloc_prop.get_subview_info(); const int idim = info.dim_idx; - const int k = info.slice_idx; - - // So far we can only subview at first or second dimension. - EKAT_REQUIRE_MSG (idim==0 || idim==1, - "Error! Subview dimension index is out of bounds.\n"); - - // Use correct subview utility - if (idim==0) { - return DstView(ekat::subview(v_np1,k)); - } else { - return DstView(ekat::subview_1(v_np1,k)); + const int k = info.slice_idx; + const int k_end = info.slice_idx_end; + + std::cout << "k_end = " << k_end << "\n"; + std::cout << "sv_alloc_prop.contiguous() = " << sv_alloc_prop.contiguous() << "\n"; + + if (k_end == -1) { + EKAT_REQUIRE_MSG( + DstRank == 1 && fl.rank() == 1, + "Error! Single-slice, strided view requires destination subview" + " rank to be equal to parent field rank.\n"); + EKAT_REQUIRE_MSG( + DstRankDynamic == 1, + "Error! Strided view not allowed with compile-time dimensions.\n"); + + // Take 2 dimensional view with normal LayoutRight + auto v_np1 = f.get_ND_view(); + + // As of now, we can only single-slice subview at first or second dimension. + EKAT_REQUIRE_MSG(idim == 0 || idim == 1, + "Error! Subview dimension index is out of bounds.\n"); + + // Use correct subview utility + if (idim == 0) { + return DstView(ekat::subview(v_np1, k)); + } else { + return DstView(ekat::subview_1(v_np1, k)); + } + // k_end has been set, so we're multi-slicing, and not contiguous, so we + } else if (k_end > 0) { + EKAT_REQUIRE_MSG(DstRank == fl.rank(), + "Error! Destination view rank must be equal to parent " + "field's rank for multi-sliced subview.\n"); + auto v_fullsize = f.get_ND_view(); + + return DstView(ekat::subview( + v_fullsize, Kokkos::make_pair(k, k_end), idim)); } } - // Not a subfield, so stride=1, and we can create the strided view from the LayoutRight 1d view. - return DstView(get_ND_view()); + // Not a subfield, so stride=1, and we can create the strided view from the + // LayoutRight 1d view. + return DstView(get_ND_view()); } -// NOTE: multi-slicing a view is only supported for strided view return type -// TODO: N doesn't actually need to be passed in, since the subview must have -// the same rank as its parent -template -auto Field::get_multi_sliced_view() const - -> get_strided_view_type, HD> { - // The destination view type on correct mem space - using DstView = get_strided_view_type, HD>; - // The dst value types - using DstValueType = typename DstView::traits::value_type; - - // Get src details - const auto& alloc_prop = m_header->get_alloc_properties(); - const auto& fl = m_header->get_identifier().get_layout(); - - // Checks - EKAT_REQUIRE_MSG(N == fl.rank(), "Error! Input Rank must be equal to parent " - "view's rank for multi-sliced subview.\n"); - EKAT_REQUIRE_MSG( - is_allocated(), - "Error! Cannot extract a field's view before allocation happens.\n"); - EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, - "Error! Cannot get a view to non-const data if the field is " - "read-only.\n"); - EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), - "Error! Source field allocation is not compatible with the " - "requested value type.\n"); - - // get parent header and check if this field is a subview of another field - const auto parent = m_header->get_parent().lock(); - EKAT_REQUIRE_MSG( - parent != nullptr, - "Error! Multi-sliced subview is unavailable for non-subfields.\n"); - Field f; - // create new field with the header/data from the parent - f.m_header = parent; - f.m_data = m_data; - - auto v_fullsize = f.get_ND_view(); - - // Now we can subview v_np1 at the correct slice - const auto& info = m_header->get_alloc_properties().get_subview_info(); - const int idim = info.dim_idx; - const int k_beg = info.slice_idx; - const int k_end = info.slice_idx_end; - - return DstView(ekat::subview( - v_fullsize, Kokkos::make_pair(k_beg, k_end), idim)); -} +// // NOTE: multi-slicing a view is only supported for strided view return type +// // TODO: N doesn't actually need to be passed in, since the subview must have +// // the same rank as its parent +// template +// auto Field::get_multi_sliced_view() const +// -> get_strided_view_type, HD> { +// // The destination view type on correct mem space +// using DstView = get_strided_view_type, HD>; +// // The dst value types +// using DstValueType = typename DstView::traits::value_type; + +// // Get src details +// const auto& alloc_prop = m_header->get_alloc_properties(); +// const auto& fl = m_header->get_identifier().get_layout(); + +// // Checks +// EKAT_REQUIRE_MSG(N == fl.rank(), "Error! Input Rank must be equal to parent " +// "view's rank for multi-sliced subview.\n"); +// EKAT_REQUIRE_MSG( +// is_allocated(), +// "Error! Cannot extract a field's view before allocation happens.\n"); +// EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, +// "Error! Cannot get a view to non-const data if the field is " +// "read-only.\n"); +// EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), +// "Error! Source field allocation is not compatible with the " +// "requested value type.\n"); + +// // get parent header and check if this field is a subview of another field +// const auto parent = m_header->get_parent().lock(); +// EKAT_REQUIRE_MSG( +// parent != nullptr, +// "Error! Multi-sliced subview is unavailable for non-subfields.\n"); +// Field f; +// // create new field with the header/data from the parent +// f.m_header = parent; +// f.m_data = m_data; + +// auto v_fullsize = f.get_ND_view(); + +// // Now we can subview v_np1 at the correct slice +// const auto& info = m_header->get_alloc_properties().get_subview_info(); +// const int idim = info.dim_idx; +// const int k_beg = info.slice_idx; +// const int k_end = info.slice_idx_end; + +// return DstView(ekat::subview( +// v_fullsize, Kokkos::make_pair(k_beg, k_end), idim)); +// } template void Field:: diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp index b68e62b4fdd..a20a945f651 100644 --- a/components/eamxx/src/share/tests/subfield_tests.cpp +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -47,32 +47,32 @@ TEST_CASE("field", "") { } SECTION("multi-sliced subfield") { - SECTION("1D multi-slice") { - // ============== - /* Rank-1 view */ - // ============== - std::vector t1 = {COL}; - std::vector d1 = {11}; - FieldIdentifier fid1("1d", {t1, d1}, m / s, "some_grid"); + // SECTION("1D multi-slice") { + // // ============== + // /* Rank-1 view */ + // // ============== + // std::vector t1 = {COL}; + // std::vector d1 = {11}; + // FieldIdentifier fid1("1d", {t1, d1}, m / s, "some_grid"); - Field f1(fid1); - f1.allocate_view(); - randomize(f1, engine, pdf); + // Field f1(fid1); + // f1.allocate_view(); + // randomize(f1, engine, pdf); - const int idim = {0}; - const int sl_beg = {3}; - const int sl_end = {7}; + // const int idim = {0}; + // const int sl_beg = {3}; + // const int sl_end = {7}; - auto v1d_h = f1.get_view(); - auto sf = f1.subfield(idim, sl_beg, sl_end); - auto sv_h = sf.get_multi_sliced_view(); + // auto v1d_h = f1.get_view(); + // auto sf = f1.subfield(idim, sl_beg, sl_end); + // auto sv_h = sf.get_strided_view(); - REQUIRE(sv_h.extent_int(idim) == (sl_end - sl_beg)); + // REQUIRE(sv_h.extent_int(idim) == (sl_end - sl_beg)); - for (int i = sl_beg; i < sl_end; i++) { - REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); - } - } + // for (int i = sl_beg; i < sl_end; i++) { + // REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); + // } + // } SECTION("2D multi-slice") { std::vector t2 = {COL, CMP}; @@ -92,7 +92,7 @@ TEST_CASE("field", "") { for (int ens = 0; ens < 2; ens++) { auto sf = f2.subfield(idim[ens], sl_beg[ens], sl_end[ens]); - auto sv_h = sf.get_multi_sliced_view(); + auto sv_h = sf.get_strided_view(); i1 = (ens == 0) ? sl_beg[0] : 0; i2 = (ens == 0) ? sl_end[0] : d2[0]; j1 = (ens == 1) ? sl_beg[1] : 0; @@ -125,7 +125,7 @@ TEST_CASE("field", "") { for (int ens = 0; ens < 3; ens++) { auto sf = f3.subfield(idim[ens], sl_beg[ens], sl_end[ens]); - auto sv_h = sf.get_multi_sliced_view(); + auto sv_h = sf.get_strided_view(); i1 = (ens == 0) ? sl_beg[0] : 0; i2 = (ens == 0) ? sl_end[0] : d3[0]; j1 = (ens == 1) ? sl_beg[1] : 0; @@ -145,7 +145,7 @@ TEST_CASE("field", "") { // ====================== auto cmp3 = f3.get_components(sl_beg[1], sl_end[1]); - auto svc_h = cmp3.get_multi_sliced_view(); + auto svc_h = cmp3.get_strided_view(); for (int i = 0; i < d3[0]; i++) { for (int j = sl_beg[1]; j < sl_end[1]; j++) { for (int k = 0; k < d3[2]; k++) { @@ -175,7 +175,7 @@ TEST_CASE("field", "") { for (int ens = 0; ens < 4; ens++) { auto sf = f4.subfield(idim[ens], sl_beg[ens], sl_end[ens]); - auto sv_h = sf.get_multi_sliced_view(); + auto sv_h = sf.get_strided_view(); i1 = (ens == 0) ? sl_beg[0] : 0; i2 = (ens == 0) ? sl_end[0] : d4[0]; j1 = (ens == 1) ? sl_beg[1] : 0; @@ -217,7 +217,7 @@ TEST_CASE("field", "") { for (int ens = 0; ens < 5; ens++) { auto sf = f5.subfield(idim[ens], sl_beg[ens], sl_end[ens]); - auto sv_h = sf.get_multi_sliced_view(); + auto sv_h = sf.get_strided_view(); i1 = (ens == 0) ? sl_beg[0] : 0; i2 = (ens == 0) ? sl_end[0] : d5[0]; j1 = (ens == 1) ? sl_beg[1] : 0; @@ -263,7 +263,7 @@ TEST_CASE("field", "") { for (int ens = 0; ens < 6; ens++) { auto sf = f6.subfield(idim[ens], sl_beg[ens], sl_end[ens]); - auto sv_h = sf.get_multi_sliced_view(); + auto sv_h = sf.get_strided_view(); i1 = (ens == 0) ? sl_beg[0] : 0; i2 = (ens == 0) ? sl_end[0] : d6[0]; j1 = (ens == 1) ? sl_beg[1] : 0; From 950e19b7ff7adb34a1e774cdfb400a7185a7600f Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Thu, 25 Apr 2024 13:18:38 -0600 Subject: [PATCH 10/16] unify multi-slice subviews with single-slice under get_strided_view() --- components/eamxx/src/share/field/field.hpp | 25 +++-- .../eamxx/src/share/field/field_impl.hpp | 96 +++++++------------ .../src/share/field/field_utils_impl.hpp | 33 ++++--- .../eamxx/src/share/tests/subfield_tests.cpp | 55 +++++------ externals/ekat | 2 +- 5 files changed, 101 insertions(+), 110 deletions(-) diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index f4ad8199726..5b3935fba79 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -318,25 +318,38 @@ class Field { if_t<(N<=2), get_view_type,HD>> get_subview_1 (const get_view_type,HD>&, const int) const { - EKAT_ERROR_MSG ("Error! Cannot subview a rank2 view along the second dimension without losing LayoutRight.\n"); + EKAT_ERROR_MSG ("Error! Cannot subview a rank2 view along the second " + "dimension without losing LayoutRight.\n"); } template auto get_ND_view () const - -> if_t,HD>>; + -> if_t<((N < MaxRank) and (N > 0)), get_view_type,HD>>; template auto get_ND_view () const - -> if_t<(N,HD>>; + -> if_t,HD>>; + + // NOTE: DO NOT USE--this circumvents compile-time issues with + // subview slicing in get_strided_view() + template + auto get_ND_view () const + -> if_t,HD>>; + + // NOTE: DO NOT USE--it only returns an error and is here to protect + // against compiler errors related to sliced subviews in get_strided_view() + template + auto get_ND_view () const + -> if_t<(N >= MaxRank + 1), get_view_type,HD>>; // Metadata (name, rank, dims, customer/providers, time stamp, ...) - std::shared_ptr m_header; + std::shared_ptr m_header; // Actual data. - dual_view_t m_data; + dual_view_t m_data; // Whether this field is read-only - bool m_is_read_only = false; + bool m_is_read_only = false; }; // We use this to find a Field in a std container. diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index c2fdce5d5cf..f93a49c809e 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -90,7 +90,8 @@ auto Field::get_view () const using DstView = get_view_type; // The dst value types using DstValueType = typename DstView::traits::value_type; - // The ViewDimension object from the Dst View (used to check validity of possible compile-time extents) + // The ViewDimension object from the Dst View (used to check validity of + // possible compile-time extents) using dims_type = typename DstView::traits::dimension; // We only allow to reshape to a view of the correct rank constexpr int DstRank = DstView::rank; @@ -185,9 +186,7 @@ get_strided_view_type { const int k = info.slice_idx; const int k_end = info.slice_idx_end; - std::cout << "k_end = " << k_end << "\n"; - std::cout << "sv_alloc_prop.contiguous() = " << sv_alloc_prop.contiguous() << "\n"; - + // k_end has not been set by a multi-slice subfield function if (k_end == -1) { EKAT_REQUIRE_MSG( DstRank == 1 && fl.rank() == 1, @@ -197,8 +196,9 @@ get_strided_view_type { DstRankDynamic == 1, "Error! Strided view not allowed with compile-time dimensions.\n"); - // Take 2 dimensional view with normal LayoutRight - auto v_np1 = f.get_ND_view(); + // Take an (n + 1)-dimensional == DstRank (== 2D, in practice) view + // with normal LayoutRight + auto v_np1 = f.get_ND_view(); // As of now, we can only single-slice subview at first or second dimension. EKAT_REQUIRE_MSG(idim == 0 || idim == 1, @@ -210,7 +210,7 @@ get_strided_view_type { } else { return DstView(ekat::subview_1(v_np1, k)); } - // k_end has been set, so we're multi-slicing, and not contiguous, so we + // k_end has been set, so we're multi-slicing } else if (k_end > 0) { EKAT_REQUIRE_MSG(DstRank == fl.rank(), "Error! Destination view rank must be equal to parent " @@ -221,62 +221,11 @@ get_strided_view_type { v_fullsize, Kokkos::make_pair(k, k_end), idim)); } } - - // Not a subfield, so stride=1, and we can create the strided view from the + // Not a subfield, so stride == 1, and we can create the strided view from the // LayoutRight 1d view. - return DstView(get_ND_view()); + return DstView(get_ND_view()); } -// // NOTE: multi-slicing a view is only supported for strided view return type -// // TODO: N doesn't actually need to be passed in, since the subview must have -// // the same rank as its parent -// template -// auto Field::get_multi_sliced_view() const -// -> get_strided_view_type, HD> { -// // The destination view type on correct mem space -// using DstView = get_strided_view_type, HD>; -// // The dst value types -// using DstValueType = typename DstView::traits::value_type; - -// // Get src details -// const auto& alloc_prop = m_header->get_alloc_properties(); -// const auto& fl = m_header->get_identifier().get_layout(); - -// // Checks -// EKAT_REQUIRE_MSG(N == fl.rank(), "Error! Input Rank must be equal to parent " -// "view's rank for multi-sliced subview.\n"); -// EKAT_REQUIRE_MSG( -// is_allocated(), -// "Error! Cannot extract a field's view before allocation happens.\n"); -// EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, -// "Error! Cannot get a view to non-const data if the field is " -// "read-only.\n"); -// EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), -// "Error! Source field allocation is not compatible with the " -// "requested value type.\n"); - -// // get parent header and check if this field is a subview of another field -// const auto parent = m_header->get_parent().lock(); -// EKAT_REQUIRE_MSG( -// parent != nullptr, -// "Error! Multi-sliced subview is unavailable for non-subfields.\n"); -// Field f; -// // create new field with the header/data from the parent -// f.m_header = parent; -// f.m_data = m_data; - -// auto v_fullsize = f.get_ND_view(); - -// // Now we can subview v_np1 at the correct slice -// const auto& info = m_header->get_alloc_properties().get_subview_info(); -// const int idim = info.dim_idx; -// const int k_beg = info.slice_idx; -// const int k_end = info.slice_idx_end; - -// return DstView(ekat::subview( -// v_fullsize, Kokkos::make_pair(k_beg, k_end), idim)); -// } - template void Field:: deep_copy (const Field& src) { @@ -343,7 +292,7 @@ deep_copy_impl (const Field& src) { const auto& layout_src = src.get_header().get_identifier().get_layout(); EKAT_REQUIRE_MSG(layout==layout_src, "ERROR: Unable to copy field " + src.get_header().get_identifier().name() + - " to field " + get_header().get_identifier().name() + ". Layouts don't match."); + " to field " + get_header().get_identifier().name() + ". Layouts don't match.\n"); const auto rank = layout.rank(); // For rank 0 view, we only need to copy a single value and return @@ -778,7 +727,7 @@ update_impl (const Field& x, const ST alpha, const ST beta, const ST fill_val) template auto Field::get_ND_view () const -> - if_t<(N,HD>> + if_t<((N < MaxRank) and (N > 0)),get_view_type,HD>> { const auto& fl = m_header->get_identifier().get_layout(); EKAT_REQUIRE_MSG (N==1 || N==fl.rank(), @@ -863,6 +812,29 @@ auto Field::get_ND_view () const -> return ret_type (ptr,kl); } +// NOTE: DO NOT USE--this circumvents compile-time issues with +// subview slicing in get_strided_view() +template +auto Field::get_ND_view () const -> + if_t,HD>> +{ + EKAT_ERROR_MSG("Error! Cannot call get_ND_view for rank == 0 (scalar).\n" + "This should never be called at run time.\n" + "Please contact developer if this functionality is required\n"); +} + +// NOTE: DO NOT USE--this circumvents compile-time issues with +// subview slicing in get_strided_view() +template +auto Field::get_ND_view () const -> + if_t<(N >= MaxRank + 1),get_view_type,HD>> +{ + EKAT_ERROR_MSG("Error! Cannot call get_ND_view for rank greater than " + "MaxRank = 6.\n" + "This should never be called at run time.\n" + "Please contact developer if this functionality is required\n"); +} + } // namespace scream #endif // SCREAM_FIELD_IMPL_HPP diff --git a/components/eamxx/src/share/field/field_utils_impl.hpp b/components/eamxx/src/share/field/field_utils_impl.hpp index 222ff1bb49e..d6ed0b13d44 100644 --- a/components/eamxx/src/share/field/field_utils_impl.hpp +++ b/components/eamxx/src/share/field/field_utils_impl.hpp @@ -29,6 +29,8 @@ bool views_are_equal(const Field& f1, const Field& f2, const ekat::Comm* comm) f2.sync_to_host(); // Reshape based on the rank, then loop over all entries. + // NOTE: because views_are_equal() is only used for testing, we generalize by + // always calling get_strided_view(), even if it could be contiguous bool same_locally = true; const auto& dims = l1.dims(); switch (l1.rank()) { @@ -46,8 +48,8 @@ bool views_are_equal(const Field& f1, const Field& f2, const ekat::Comm* comm) break; case 2: { - auto v1 = f1.template get_view(); - auto v2 = f2.template get_view(); + auto v1 = f1.template get_strided_view(); + auto v2 = f2.template get_strided_view(); for (int i=0; same_locally && i(); - auto v2 = f2.template get_view(); + auto v1 = f1.template get_strided_view(); + auto v2 = f2.template get_strided_view(); for (int i=0; same_locally && i(); - auto v2 = f2.template get_view(); + auto v1 = f1.template get_strided_view(); + auto v2 = f2.template get_strided_view(); for (int i=0; same_locally && i(); - auto v2 = f2.template get_view(); + auto v1 = f1.template get_strided_view(); + auto v2 = f2.template get_strided_view(); for (int i=0; same_locally && i(); - auto v2 = f2.template get_view(); + auto v1 = f1.template get_strided_view(); + auto v2 = f2.template get_strided_view(); for (int i=0; same_locally && i()() << ", \n"; + out << "\n " << f.get_strided_view()() << ", \n"; break; } case 1: @@ -769,7 +774,7 @@ void print_field_hyperslab (const Field& f, case 2: { dims_str[dims_left[1]] = ":"; - auto v = f.get_view(); + auto v = f.get_strided_view(); for (int i=0; i(); + auto v = f.get_strided_view(); for (int i=0; i(); + auto v = f.get_strided_view(); for (int i=0; i t1 = {COL}; - // std::vector d1 = {11}; - // FieldIdentifier fid1("1d", {t1, d1}, m / s, "some_grid"); - - // Field f1(fid1); - // f1.allocate_view(); - // randomize(f1, engine, pdf); - - // const int idim = {0}; - // const int sl_beg = {3}; - // const int sl_end = {7}; - - // auto v1d_h = f1.get_view(); - // auto sf = f1.subfield(idim, sl_beg, sl_end); - // auto sv_h = sf.get_strided_view(); - - // REQUIRE(sv_h.extent_int(idim) == (sl_end - sl_beg)); - - // for (int i = sl_beg; i < sl_end; i++) { - // REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); - // } - // } + SECTION("1D multi-slice") { + // ============== + /* Rank-1 view */ + // ============== + std::vector t1 = {COL}; + std::vector d1 = {11}; + FieldIdentifier fid1("1d", {t1, d1}, m / s, "some_grid"); + + Field f1(fid1); + f1.allocate_view(); + randomize(f1, engine, pdf); + + const int idim = {0}; + const int sl_beg = {3}; + const int sl_end = {7}; + + auto v1d_h = f1.get_view(); + auto sf = f1.subfield(idim, sl_beg, sl_end); + REQUIRE(sf.get_header().get_alloc_properties().contiguous() == false); + auto sv_h = sf.get_strided_view(); + REQUIRE(sv_h.extent_int(idim) == (sl_end - sl_beg)); + + for (int i = sl_beg; i < sl_end; i++) { + REQUIRE(v1d_h(i) == sv_h(i - sl_beg)); + } + } SECTION("2D multi-slice") { std::vector t2 = {COL, CMP}; @@ -145,7 +145,7 @@ TEST_CASE("field", "") { // ====================== auto cmp3 = f3.get_components(sl_beg[1], sl_end[1]); - auto svc_h = cmp3.get_strided_view(); + auto svc_h = cmp3.get_strided_view(); for (int i = 0; i < d3[0]; i++) { for (int j = sl_beg[1]; j < sl_end[1]; j++) { for (int k = 0; k < d3[2]; k++) { @@ -263,6 +263,7 @@ TEST_CASE("field", "") { for (int ens = 0; ens < 6; ens++) { auto sf = f6.subfield(idim[ens], sl_beg[ens], sl_end[ens]); + REQUIRE(sf.get_header().get_alloc_properties().contiguous() == false); auto sv_h = sf.get_strided_view(); i1 = (ens == 0) ? sl_beg[0] : 0; i2 = (ens == 0) ? sl_end[0] : d6[0]; diff --git a/externals/ekat b/externals/ekat index eed9ddbefe6..bb3ebcdfd5f 160000 --- a/externals/ekat +++ b/externals/ekat @@ -1 +1 @@ -Subproject commit eed9ddbefe685ea084853fa9c9262c861a5a7f39 +Subproject commit bb3ebcdfd5f600c8e75420eaed167d9b268b100b From 51dae5b253c148731cc64d8a256e6f5bd77cdb0e Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Fri, 26 Apr 2024 00:40:01 -0600 Subject: [PATCH 11/16] fix other tests in share/ --- components/eamxx/src/share/field/field.hpp | 8 +-- .../eamxx/src/share/field/field_impl.hpp | 21 +----- .../src/share/field/field_utils_impl.hpp | 64 +++++++++---------- .../share/property_checks/field_nan_check.cpp | 14 ++-- .../eamxx/src/share/tests/field_utils.cpp | 26 ++++---- .../eamxx/src/share/tests/property_checks.cpp | 12 ++-- externals/ekat | 2 +- 7 files changed, 61 insertions(+), 86 deletions(-) diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index 5b3935fba79..a3ed9076a2a 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -324,18 +324,12 @@ class Field { template auto get_ND_view () const - -> if_t<((N < MaxRank) and (N > 0)), get_view_type,HD>>; + -> if_t<(N < MaxRank), get_view_type,HD>>; template auto get_ND_view () const -> if_t,HD>>; - // NOTE: DO NOT USE--this circumvents compile-time issues with - // subview slicing in get_strided_view() - template - auto get_ND_view () const - -> if_t,HD>>; - // NOTE: DO NOT USE--it only returns an error and is here to protect // against compiler errors related to sliced subviews in get_strided_view() template diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index f93a49c809e..fcfae128ad9 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -188,14 +188,6 @@ get_strided_view_type { // k_end has not been set by a multi-slice subfield function if (k_end == -1) { - EKAT_REQUIRE_MSG( - DstRank == 1 && fl.rank() == 1, - "Error! Single-slice, strided view requires destination subview" - " rank to be equal to parent field rank.\n"); - EKAT_REQUIRE_MSG( - DstRankDynamic == 1, - "Error! Strided view not allowed with compile-time dimensions.\n"); - // Take an (n + 1)-dimensional == DstRank (== 2D, in practice) view // with normal LayoutRight auto v_np1 = f.get_ND_view(); @@ -727,7 +719,7 @@ update_impl (const Field& x, const ST alpha, const ST beta, const ST fill_val) template auto Field::get_ND_view () const -> - if_t<((N < MaxRank) and (N > 0)),get_view_type,HD>> + if_t<(N < MaxRank), get_view_type,HD>> { const auto& fl = m_header->get_identifier().get_layout(); EKAT_REQUIRE_MSG (N==1 || N==fl.rank(), @@ -812,17 +804,6 @@ auto Field::get_ND_view () const -> return ret_type (ptr,kl); } -// NOTE: DO NOT USE--this circumvents compile-time issues with -// subview slicing in get_strided_view() -template -auto Field::get_ND_view () const -> - if_t,HD>> -{ - EKAT_ERROR_MSG("Error! Cannot call get_ND_view for rank == 0 (scalar).\n" - "This should never be called at run time.\n" - "Please contact developer if this functionality is required\n"); -} - // NOTE: DO NOT USE--this circumvents compile-time issues with // subview slicing in get_strided_view() template diff --git a/components/eamxx/src/share/field/field_utils_impl.hpp b/components/eamxx/src/share/field/field_utils_impl.hpp index d6ed0b13d44..e12f6debbc1 100644 --- a/components/eamxx/src/share/field/field_utils_impl.hpp +++ b/components/eamxx/src/share/field/field_utils_impl.hpp @@ -141,13 +141,13 @@ void randomize (const Field& f, Engine& engine, PDF&& pdf) switch (fl.rank()) { case 0: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); v() = pdf(engine); } break; case 1: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + const auto gids = dof_gids.get_strided_view(); // Create a field to store perturbation values with layout // the same as f, but stripped of column and level dimension. @@ -296,7 +296,7 @@ ST frobenius_norm(const Field& f, const ekat::Comm* comm) switch (fl.rank()) { case 1: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); for (int i=0; i(); + auto v = f.template get_strided_view(); Kokkos::parallel_reduce(size, KOKKOS_LAMBDA(int i, int& result) { if (ekat::is_invalid(v(i))) { result = i; @@ -57,7 +57,7 @@ PropertyCheck::ResultAndMsg FieldNaNCheck::check_impl() const { break; case 2: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); Kokkos::parallel_reduce(size, KOKKOS_LAMBDA(int idx, int& result) { int i,j; unflatten_idx(idx,extents,i,j); @@ -69,7 +69,7 @@ PropertyCheck::ResultAndMsg FieldNaNCheck::check_impl() const { break; case 3: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); Kokkos::parallel_reduce(size, KOKKOS_LAMBDA(int idx, int& result) { int i,j,k; unflatten_idx(idx,extents,i,j,k); @@ -81,7 +81,7 @@ PropertyCheck::ResultAndMsg FieldNaNCheck::check_impl() const { break; case 4: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); Kokkos::parallel_reduce(size, KOKKOS_LAMBDA(int idx, int& result) { int i,j,k,l; unflatten_idx(idx,extents,i,j,k,l); @@ -93,7 +93,7 @@ PropertyCheck::ResultAndMsg FieldNaNCheck::check_impl() const { break; case 5: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); Kokkos::parallel_reduce(size, KOKKOS_LAMBDA(int idx, int& result) { int i,j,k,l,m; unflatten_idx(idx,extents,i,j,k,l,m); @@ -105,7 +105,7 @@ PropertyCheck::ResultAndMsg FieldNaNCheck::check_impl() const { break; case 6: { - auto v = f.template get_view(); + auto v = f.template get_strided_view(); Kokkos::parallel_reduce(size, KOKKOS_LAMBDA(int idx, int& result) { int i,j,k,l,m,n; unflatten_idx(idx,extents,i,j,k,l,m,n); @@ -136,7 +136,7 @@ PropertyCheck::ResultAndMsg FieldNaNCheck::check_impl() const { if (m_grid) { // We are storing grid info, and the field is over columns. Get col id and coords. col_lid = indices[0]; - auto gids = m_grid->get_dofs_gids().get_view(); + auto gids = m_grid->get_dofs_gids().get_strided_view(); res_and_msg.msg += " - indices (w/ global column index): (" + std::to_string(gids(col_lid)); for (size_t i=1; i(); - auto v2 = f2.get_view(); + auto v1 = f1.get_strided_view(); + auto v2 = f2.get_strided_view(); auto dim0 = fid.get_layout().dim(0); auto dim1 = fid.get_layout().dim(1); auto am_i_root = comm.am_i_root(); @@ -80,7 +80,7 @@ TEST_CASE("utils") { SECTION ("sum") { - auto v1 = f1.get_view(); + auto v1 = f1.get_strided_view(); auto dim0 = fid.get_layout().dim(0); auto dim1 = fid.get_layout().dim(1); auto lsize = fid.get_layout().size(); @@ -103,7 +103,7 @@ TEST_CASE("utils") { SECTION ("frobenius") { - auto v1 = f1.get_view(); + auto v1 = f1.get_strided_view(); auto dim0 = fid.get_layout().dim(0); auto dim1 = fid.get_layout().dim(1); auto lsize = fid.get_layout().size(); @@ -130,7 +130,7 @@ TEST_CASE("utils") { SECTION ("max") { - auto v1 = f1.get_view(); + auto v1 = f1.get_strided_view(); auto dim0 = fid.get_layout().dim(0); auto dim1 = fid.get_layout().dim(1); auto lsize = fid.get_layout().size(); @@ -153,7 +153,7 @@ TEST_CASE("utils") { SECTION ("min") { - auto v1 = f1.get_view(); + auto v1 = f1.get_strided_view(); auto dim0 = fid.get_layout().dim(0); auto dim1 = fid.get_layout().dim(1); auto lsize = fid.get_layout().size(); @@ -221,10 +221,10 @@ TEST_CASE("utils") { // Sync to host for checks f1.sync_to_host(), f2a.sync_to_host(), f2b.sync_to_host(), f3.sync_to_host(); - const auto v1 = f1.get_view (); - const auto v2a = f2a.get_view(); - const auto v2b = f2b.get_view(); - const auto v3 = f3.get_view (); + const auto v1 = f1.get_strided_view (); + const auto v2a = f2a.get_strided_view(); + const auto v2b = f2b.get_strided_view(); + const auto v3 = f3.get_strided_view (); // Check that all field values are 1 for all but last 3 levels and between [2,3] otherwise. auto check_level = [&] (const int ilev, const Real val) { @@ -248,8 +248,8 @@ TEST_CASE("utils") { perturb(f3_alt, engine, pdf, base_seed_alt, mask_lambda, gids); f1_alt.sync_to_host(), f3_alt.sync_to_host(); - const auto v1_alt = f1_alt.get_view(); - const auto v3_alt = f3_alt.get_view(); + const auto v1_alt = f1_alt.get_strided_view(); + const auto v3_alt = f3_alt.get_strided_view(); auto check_diff = [&] (const int ilev, const Real val1, const Real val2) { if (ilev < nlevs-3) REQUIRE(val1==val2); @@ -318,7 +318,7 @@ TEST_CASE ("print_field_hyperslab") { f.allocate_view(); randomize (f,engine,pdf); - auto v = f.get_view(); + auto v = f.get_strided_view(); SECTION ("slice_0") { std::vector loc_tags = {EL,CMP}; diff --git a/components/eamxx/src/share/tests/property_checks.cpp b/components/eamxx/src/share/tests/property_checks.cpp index 1cb5e64369e..e0da898a841 100644 --- a/components/eamxx/src/share/tests/property_checks.cpp +++ b/components/eamxx/src/share/tests/property_checks.cpp @@ -69,10 +69,10 @@ TEST_CASE("property_checks", "") { const auto units = ekat::units::Units::nondimensional(); const auto& lat = grid->create_geometry_data("lat",layout,units); const auto& lon = grid->create_geometry_data("lon",layout,units); - auto lat_h = lat.get_view(); - auto lon_h = lon.get_view(); + auto lat_h = lat.get_strided_view(); + auto lon_h = lon.get_strided_view(); auto dofs = grid->get_dofs_gids(); - auto dofs_h = dofs.get_view(); + auto dofs_h = dofs.get_strided_view(); for (int i=0; iget_num_local_dofs(); ++i) { lat_h(i) = i; lon_h(i) = -i; @@ -110,12 +110,12 @@ TEST_CASE("property_checks", "") { REQUIRE(res_and_msg.result==CheckResult::Pass); // Assign a NaN value to the field, make sure it fails the check, - auto f_view = f.get_view(); + auto f_view = f.get_strided_view(); f_view(1,2,3) = std::numeric_limits::quiet_NaN(); f.sync_to_dev(); res_and_msg = nan_check->check(); REQUIRE(res_and_msg.result==CheckResult::Fail); - + std::string expected_msg = "FieldNaNCheck failed.\n" " - field id: " + fid.get_id_string() + "\n" @@ -149,7 +149,7 @@ TEST_CASE("property_checks", "") { // Assign out-of-bounds values to the field, make sure it fails the check, // and then repair the field so it passes. f.deep_copy(0.5); - auto f_view = f.get_view(); + auto f_view = f.get_strided_view(); f_view(1,2,3) = 2.0; f_view(0,1,2) = 0.0; f.sync_to_dev(); diff --git a/externals/ekat b/externals/ekat index bb3ebcdfd5f..c5c05e62664 160000 --- a/externals/ekat +++ b/externals/ekat @@ -1 +1 @@ -Subproject commit bb3ebcdfd5f600c8e75420eaed167d9b268b100b +Subproject commit c5c05e626648fa62fff8d99c51457986acdbc0c4 From acf6214ddeb7b84767c720fb63048721c6ffcc02 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Tue, 21 May 2024 00:17:47 -0600 Subject: [PATCH 12/16] small changes to fix rebase --- components/eamxx/src/share/field/field.cpp | 7 +++---- components/eamxx/src/share/field/field.hpp | 6 +++++- components/eamxx/src/share/field/field_alloc_prop.cpp | 4 ++-- components/eamxx/src/share/field/field_impl.hpp | 9 ++++++++- 4 files changed, 18 insertions(+), 8 deletions(-) diff --git a/components/eamxx/src/share/field/field.cpp b/components/eamxx/src/share/field/field.cpp index f56b3b83d93..e8a3ab4a058 100644 --- a/components/eamxx/src/share/field/field.cpp +++ b/components/eamxx/src/share/field/field.cpp @@ -128,8 +128,8 @@ Field Field::subfield(const std::string& sf_name, is_allocated(), "Error! Input field must be allocated in order to subview it.\n"); - auto sf_layout = - lt.clone_with_different_extent(idim, index_end - index_beg); + auto sf_layout = lt.clone(); + sf_layout.reset_dim(idim, index_end - index_beg); // Create identifier for subfield FieldIdentifier sf_id(sf_name, sf_layout, sf_units, id.get_grid_name()); @@ -178,8 +178,7 @@ Field Field::get_components(const int i1, const int i2) { EKAT_REQUIRE_MSG(layout.is_vector_layout(), "Error! 'get_component' available only for vector fields.\n" " Layout of '" + - fname + "': " + e2str(get_layout_type(layout.tags())) + - "\n"); + fname + "': " + e2str(layout.type()) + "\n"); const int idim = layout.get_vector_component_idx(); EKAT_REQUIRE_MSG(i1 >= 0 && i2 < layout.dim(idim), diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index a3ed9076a2a..c0a84b10130 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -324,7 +324,11 @@ class Field { template auto get_ND_view () const - -> if_t<(N < MaxRank), get_view_type,HD>>; + -> if_t, HD>>; + + template + auto get_ND_view () const + -> if_t<(N > 0) and (N < MaxRank), get_view_type,HD>>; template auto get_ND_view () const diff --git a/components/eamxx/src/share/field/field_alloc_prop.cpp b/components/eamxx/src/share/field/field_alloc_prop.cpp index 635e62bae19..48796a826bd 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.cpp +++ b/components/eamxx/src/share/field/field_alloc_prop.cpp @@ -95,8 +95,8 @@ FieldAllocProp FieldAllocProp::subview(const int idim, FieldAllocProp props(m_scalar_type_size); props.m_committed = true; props.m_scalar_type_size = m_scalar_type_size; - props.m_layout = - m_layout.clone_with_different_extent(idim, k_end - k_beg); + props.m_layout = m_layout.clone(); + props.m_layout.reset_dim(idim, k_end - k_beg); // Output is contiguous if either // - this->m_contiguous=true AND idim==0 diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index fcfae128ad9..97cb2c6b36a 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -717,9 +717,16 @@ update_impl (const Field& x, const ST alpha, const ST beta, const ST fill_val) Kokkos::fence(); } +template +auto get_ND_view () const + -> if_t,HD>> +{ + EKAT_ERROR_MSG("Error! Cannot take an ND view of 0 dimensional field.\n"); +} + template auto Field::get_ND_view () const -> - if_t<(N < MaxRank), get_view_type,HD>> + if_t<(N > 0) and (N < MaxRank), get_view_type, HD>> { const auto& fl = m_header->get_identifier().get_layout(); EKAT_REQUIRE_MSG (N==1 || N==fl.rank(), From ca0b2347ff53d9b363fe2d23335cb5aaba8414e6 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Wed, 22 May 2024 17:56:03 -0600 Subject: [PATCH 13/16] discard changes on ekat and point to current commit on master --- externals/ekat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/externals/ekat b/externals/ekat index c5c05e62664..eed9ddbefe6 160000 --- a/externals/ekat +++ b/externals/ekat @@ -1 +1 @@ -Subproject commit c5c05e626648fa62fff8d99c51457986acdbc0c4 +Subproject commit eed9ddbefe685ea084853fa9c9262c861a5a7f39 From 0a18d7e123bc3185c5621110132aa77da5393af6 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Wed, 22 May 2024 18:12:48 -0600 Subject: [PATCH 14/16] update get_strided_view() approach--all appears to be working --- components/eamxx/src/share/field/field.hpp | 12 +- .../eamxx/src/share/field/field_impl.hpp | 148 +++++++++--------- 2 files changed, 74 insertions(+), 86 deletions(-) diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index c0a84b10130..858c542fd4a 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -136,12 +136,6 @@ class Field { get_strided_view_type get_strided_view () const; - // this is the same as above but for a multi-sliced view - // which is to say, also a strided view - template - auto get_multi_sliced_view () const - -> get_strided_view_type, HD>; - // These two getters are convenience function for commonly accessed metadata. // The same info can be extracted from the metadata stored in the FieldHeader DataType data_type () const { return get_header().get_identifier().data_type(); } @@ -324,11 +318,7 @@ class Field { template auto get_ND_view () const - -> if_t, HD>>; - - template - auto get_ND_view () const - -> if_t<(N > 0) and (N < MaxRank), get_view_type,HD>>; + -> if_t<(N < MaxRank), get_view_type,HD>>; template auto get_ND_view () const diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index 97cb2c6b36a..420d90492e7 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -151,70 +151,74 @@ get_strided_view_type { using DstValueType = typename DstView::traits::value_type; // We only allow to reshape to a view of the correct rank constexpr int DstRank = DstView::rank; - constexpr int DstRankDynamic = DstView::rank_dynamic; - // Get src details - const auto& alloc_prop = m_header->get_alloc_properties(); - const auto& fl = m_header->get_identifier().get_layout(); - - // Checks - EKAT_REQUIRE_MSG( - is_allocated(), - "Error! Cannot extract a field's view before allocation happens.\n"); - EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, - "Error! Cannot get a view to non-const data if the field is " - "read-only.\n"); - EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), - "Error! Source field allocation is not compatible with the " - "requested value type.\n"); - - // Check if this field is a subview of another field - const auto parent = m_header->get_parent().lock(); - if (parent != nullptr) { - // Parent field has correct layout to reinterpret the view into N+1-dim view - // So create the parent field on the fly, use it to get the N+1-dim view, - // then subview it. NOTE: we can set protected members, since f is the same - // type of this class. - Field f; - f.m_header = parent; - f.m_data = m_data; - - // get subview info to determine whether we are single- or multi-slicing - const auto& sv_alloc_prop = m_header->get_alloc_properties(); - const auto& info = sv_alloc_prop.get_subview_info(); - const int idim = info.dim_idx; - const int k = info.slice_idx; - const int k_end = info.slice_idx_end; - - // k_end has not been set by a multi-slice subfield function - if (k_end == -1) { - // Take an (n + 1)-dimensional == DstRank (== 2D, in practice) view - // with normal LayoutRight - auto v_np1 = f.get_ND_view(); - - // As of now, we can only single-slice subview at first or second dimension. - EKAT_REQUIRE_MSG(idim == 0 || idim == 1, - "Error! Subview dimension index is out of bounds.\n"); - - // Use correct subview utility - if (idim == 0) { - return DstView(ekat::subview(v_np1, k)); - } else { - return DstView(ekat::subview_1(v_np1, k)); + if constexpr (DstRank > 0) { + // Get src details + const auto& alloc_prop = m_header->get_alloc_properties(); + const auto& fl = m_header->get_identifier().get_layout(); + + // Checks + EKAT_REQUIRE_MSG( + is_allocated(), + "Error! Cannot extract a field's view before allocation happens.\n"); + EKAT_REQUIRE_MSG(not m_is_read_only || std::is_const::value, + "Error! Cannot get a view to non-const data if the field is " + "read-only.\n"); + EKAT_REQUIRE_MSG(alloc_prop.template is_compatible(), + "Error! Source field allocation is not compatible with the " + "requested value type.\n"); + + // Check if this field is a subview of another field + const auto parent = m_header->get_parent().lock(); + if (parent != nullptr) { + // Parent field has correct layout to reinterpret the view into N+1-dim view, + // for single-slice subfield, and N-dim view for multi-slice subfield. + // So create the parent field on the fly, use it to get the N+{1,0}-dim view, + // then subview it. NOTE: we can set protected members, since f is the same + // type of this class. + Field f; + f.m_header = parent; + f.m_data = m_data; + + // get subview info to determine whether we are single- or multi-slicing + const auto& sv_alloc_prop = m_header->get_alloc_properties(); + const auto& info = sv_alloc_prop.get_subview_info(); + const int idim = info.dim_idx; + const int k = info.slice_idx; + const int k_end = info.slice_idx_end; + + // k_end has not been set by a multi-slice subfield function + if (k_end == -1) { + // Take an (n + 1)-dimensional == DstRank (== 2D, in practice) view + // with normal LayoutRight + auto v_np1 = f.get_ND_view(); + + // As of now, we can only single-slice subview at first or second dimension. + EKAT_REQUIRE_MSG(idim == 0 || idim == 1, + "Error! Subview dimension index is out of bounds.\n"); + + // Use correct subview utility + if (idim == 0) { + return DstView(ekat::subview(v_np1, k)); + } else { + return DstView(ekat::subview_1(v_np1, k)); + } + // k_end has been set, so we're multi-slicing + } else if (k_end > 0) { + // rank doesn't change for multi-slice + EKAT_REQUIRE_MSG(DstRank == fl.rank(), + "Error! Destination view rank must be equal to parent " + "field's rank for multi-sliced subview.\n"); + auto v_fullsize = f.get_ND_view(); + + return DstView(ekat::subview( + v_fullsize, Kokkos::make_pair(k, k_end), idim)); } - // k_end has been set, so we're multi-slicing - } else if (k_end > 0) { - EKAT_REQUIRE_MSG(DstRank == fl.rank(), - "Error! Destination view rank must be equal to parent " - "field's rank for multi-sliced subview.\n"); - auto v_fullsize = f.get_ND_view(); - - return DstView(ekat::subview( - v_fullsize, Kokkos::make_pair(k, k_end), idim)); } } - // Not a subfield, so stride == 1, and we can create the strided view from the - // LayoutRight 1d view. + // Either not a subfield or requesting a zero-D view from a + // 0D or 1D subfield, so stride == 1, and we can create the + // strided view from the LayoutRight 1d view. return DstView(get_ND_view()); } @@ -717,16 +721,9 @@ update_impl (const Field& x, const ST alpha, const ST beta, const ST fill_val) Kokkos::fence(); } -template -auto get_ND_view () const - -> if_t,HD>> -{ - EKAT_ERROR_MSG("Error! Cannot take an ND view of 0 dimensional field.\n"); -} - -template -auto Field::get_ND_view () const -> - if_t<(N > 0) and (N < MaxRank), get_view_type, HD>> +template +auto Field::get_ND_view () const + -> if_t<(N < MaxRank), get_view_type, HD>> { const auto& fl = m_header->get_identifier().get_layout(); EKAT_REQUIRE_MSG (N==1 || N==fl.rank(), @@ -754,7 +751,8 @@ auto Field::get_ND_view () const -> "Error! Subview dimension index is out of bounds.\n"); EKAT_REQUIRE_MSG (idim==0 || N>1, - "Error! Cannot subview a rank-2 (or less) view along 2nd dimension without losing LayoutRight.\n"); + "Error! Cannot subview a rank-2 (or less) view along 2nd dimension " + "without losing LayoutRight.\n"); // Use SFINAE-ed get_subview helper function to pick correct // subview impl. If N+1<=2 and idim!=0, the code craps out in the check above. @@ -785,8 +783,8 @@ auto Field::get_ND_view () const -> } template -auto Field::get_ND_view () const -> - if_t,HD>> +auto Field::get_ND_view () const + -> if_t,HD>> { const auto& fl = m_header->get_identifier().get_layout(); EKAT_REQUIRE_MSG (N==1 || N==fl.rank(), @@ -814,8 +812,8 @@ auto Field::get_ND_view () const -> // NOTE: DO NOT USE--this circumvents compile-time issues with // subview slicing in get_strided_view() template -auto Field::get_ND_view () const -> - if_t<(N >= MaxRank + 1),get_view_type,HD>> +auto Field::get_ND_view () const + -> if_t<(N >= MaxRank + 1),get_view_type,HD>> { EKAT_ERROR_MSG("Error! Cannot call get_ND_view for rank greater than " "MaxRank = 6.\n" From d0e870838358f74acc6f99a0c3ce417d5db2ad25 Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Tue, 28 May 2024 17:05:05 -0600 Subject: [PATCH 15/16] fix Luca PR comments --- components/eamxx/src/share/field/field.cpp | 16 +++++++------- components/eamxx/src/share/field/field.hpp | 4 ++-- .../src/share/field/field_alloc_prop.cpp | 21 ++++++++----------- .../src/share/field/field_utils_impl.hpp | 8 +++++-- .../share/property_checks/field_nan_check.cpp | 2 ++ 5 files changed, 28 insertions(+), 23 deletions(-) diff --git a/components/eamxx/src/share/field/field.cpp b/components/eamxx/src/share/field/field.cpp index e8a3ab4a058..320e10a6027 100644 --- a/components/eamxx/src/share/field/field.cpp +++ b/components/eamxx/src/share/field/field.cpp @@ -115,7 +115,9 @@ subfield (const int idim, const int index, const bool dynamic) const { return subfield(m_header->get_identifier().name(),idim,index,dynamic); } -// slice at index idim, entries \in [index_beg, index_end] +// slice at index idim, extracting the N = (index_end - index_beg) entries +// written in math notation: [index_beg, index_end) +// or equivalently, subF = F(index_beg, ... , index_beg + N) Field Field::subfield(const std::string& sf_name, const ekat::units::Units& sf_units, const int idim, const int index_beg, const int index_end) const { @@ -172,7 +174,7 @@ get_component (const int i, const bool dynamic) { return subfield (fname + "_" + std::to_string(i),idim,i,dynamic); } -Field Field::get_components(const int i1, const int i2) { +Field Field::get_components(const int beg, const int end) { const auto& layout = get_header().get_identifier().get_layout(); const auto& fname = get_header().get_identifier().name(); EKAT_REQUIRE_MSG(layout.is_vector_layout(), @@ -181,15 +183,15 @@ Field Field::get_components(const int i1, const int i2) { fname + "': " + e2str(layout.type()) + "\n"); const int idim = layout.get_vector_component_idx(); - EKAT_REQUIRE_MSG(i1 >= 0 && i2 < layout.dim(idim), + EKAT_REQUIRE_MSG(beg >= 0 && end < layout.dim(idim), "Error! Component index range out of bounds [0," + std::to_string(layout.dim(idim)) + ").\n"); - EKAT_REQUIRE_MSG(i1 < i2, "Error! Invalid component indices (i1 >= i2).\n"); + EKAT_REQUIRE_MSG(beg < end, "Error! Invalid component indices (beg >= end).\n"); - // Add _$i1-$i2 to the field name, to avoid issues if the subfield is stored + // Add _$beg-$end to the field name, to avoid issues if the subfield is stored // in some structure that requires unique names (e.g., a remapper) - return subfield(fname + "_" + std::to_string(i1) + "-" + std::to_string(i2), - idim, i1, i2); + return subfield(fname + "_" + std::to_string(beg) + "-" + std::to_string(end), + idim, beg, end); } bool Field::equivalent(const Field& rhs) const diff --git a/components/eamxx/src/share/field/field.hpp b/components/eamxx/src/share/field/field.hpp index 858c542fd4a..d0966daef92 100644 --- a/components/eamxx/src/share/field/field.hpp +++ b/components/eamxx/src/share/field/field.hpp @@ -248,13 +248,13 @@ class Field { const int idim, const int index_beg, const int index_end) const; Field subfield (const std::string& sf_name, const int idim, const int index_beg, const int index_end) const; - Field subfield (const int idim, const int k_beg, const int k_end) const; + Field subfield (const int idim, const int index_beg, const int index_end) const; // If this field is a vector field, get a subfield for the ith component. // If dynamic = true, it is possible to "reset" the component index at runtime. // Note: throws if this is not a vector field. Field get_component (const int i, const bool dynamic = false); // version for slicing across multiple, contiguous indices - Field get_components (const int i1, const int i2); + Field get_components (const int beg, const int end); // Checks whether the underlying view has been already allocated. bool is_allocated () const { return m_data.d_view.data()!=nullptr; } diff --git a/components/eamxx/src/share/field/field_alloc_prop.cpp b/components/eamxx/src/share/field/field_alloc_prop.cpp index 48796a826bd..266fcb88779 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.cpp +++ b/components/eamxx/src/share/field/field_alloc_prop.cpp @@ -78,17 +78,17 @@ subview (const int idim, const int k, const bool dynamic) const { } FieldAllocProp FieldAllocProp::subview(const int idim, - const int k_beg, - const int k_end) const { + const int index_beg, + const int index_end) const { EKAT_REQUIRE_MSG( is_committed(), "Error! Subview requires alloc properties to be committed.\n"); EKAT_REQUIRE_MSG(idim < m_layout.rank(), "Error! Dimension index out of bounds.\n"); - EKAT_REQUIRE_MSG(k_beg < k_end, + EKAT_REQUIRE_MSG(index_beg < index_end, "Error! Slice indices are invalid (non-increasing).\n"); EKAT_REQUIRE_MSG( - k_beg >= 0 && k_end < m_layout.dim(idim), + index_beg >= 0 && index_end < m_layout.dim(idim), "Error! Slice index range along the idim dimension is out of bounds.\n"); // Set new layout basic stuff @@ -96,23 +96,20 @@ FieldAllocProp FieldAllocProp::subview(const int idim, props.m_committed = true; props.m_scalar_type_size = m_scalar_type_size; props.m_layout = m_layout.clone(); - props.m_layout.reset_dim(idim, k_end - k_beg); + props.m_layout.reset_dim(idim, index_end - index_beg); - // Output is contiguous if either - // - this->m_contiguous=true AND idim==0 - // - m_layout.dim(i)==1 for all i; + // below, we can't be sure the field we consider has a continuous allocation, + // so we use get_strided_view() switch (layout.rank()) { case 1: { From b41796866033260290d95f561e2741180dfe4d9b Mon Sep 17 00:00:00 2001 From: Michael J Schmidt Date: Fri, 31 May 2024 08:20:40 -0600 Subject: [PATCH 16/16] address fieldAllocProp::subview() comments, and add deep_copy test for subfield --- .../src/share/field/field_alloc_prop.cpp | 12 +- .../eamxx/src/share/field/field_impl.hpp | 118 ++++++++++++------ .../eamxx/src/share/tests/subfield_tests.cpp | 36 ++++++ 3 files changed, 120 insertions(+), 46 deletions(-) diff --git a/components/eamxx/src/share/field/field_alloc_prop.cpp b/components/eamxx/src/share/field/field_alloc_prop.cpp index 266fcb88779..be07af01445 100644 --- a/components/eamxx/src/share/field/field_alloc_prop.cpp +++ b/components/eamxx/src/share/field/field_alloc_prop.cpp @@ -67,13 +67,13 @@ subview (const int idim, const int k, const bool dynamic) const { // and there is no packing props.m_last_extent = m_layout.dim(idim); props.m_pack_size_max = 1; - props.m_alloc_size = m_alloc_size / m_last_extent; } else { // We are keeping the last dim, so same last extent and max pack size props.m_last_extent = m_last_extent; props.m_pack_size_max = m_pack_size_max; - props.m_alloc_size = m_alloc_size / m_layout.dim(idim); } + // give it an invalid value because subviews don't get allocated + props.m_alloc_size = -1; return props; } @@ -103,21 +103,19 @@ FieldAllocProp FieldAllocProp::subview(const int idim, props.m_subview_info = SubviewInfo(idim, index_beg, index_end, m_layout.dim(idim)); - // Figure out strides/packs if (idim == (m_layout.rank() - 1)) { // We're slicing the possibly padded dim, so everything else is as in the // layout, and there is no packing - // props.m_last_extent = m_layout.dim(idim); props.m_last_extent = index_end - index_beg; props.m_pack_size_max = 1; - props.m_alloc_size = m_alloc_size / m_last_extent; } else { // We are keeping the last dim, so same last extent and max pack size props.m_last_extent = m_last_extent; props.m_pack_size_max = m_pack_size_max; - props.m_alloc_size = m_alloc_size / m_layout.dim(idim); } + // give it an invalid value because subviews don't get allocated + props.m_alloc_size = -1; return props; } @@ -125,7 +123,7 @@ void FieldAllocProp::request_allocation (const int pack_size) { using ekat::ScalarTraits; EKAT_REQUIRE_MSG(!m_committed, - "Error! Cannot change allocation properties after they have been commited.\n"); + "Error! Cannot change allocation properties after they have been committed.\n"); const int vts = m_scalar_type_size*pack_size; diff --git a/components/eamxx/src/share/field/field_impl.hpp b/components/eamxx/src/share/field/field_impl.hpp index 420d90492e7..42a0a6cd61e 100644 --- a/components/eamxx/src/share/field/field_impl.hpp +++ b/components/eamxx/src/share/field/field_impl.hpp @@ -305,8 +305,6 @@ deep_copy_impl (const Field& src) { // different pack sizes. auto src_alloc_props = src.get_header().get_alloc_properties(); auto tgt_alloc_props = get_header().get_alloc_properties(); - auto src_alloc_size = src_alloc_props.get_alloc_size(); - auto tgt_alloc_size = tgt_alloc_props.get_alloc_size(); // If a manual parallel_for is required (b/c of alloc sizes difference), // we need to create extents (rather than just using the one in layout), @@ -325,33 +323,32 @@ deep_copy_impl (const Field& src) { if (src_alloc_props.contiguous() and tgt_alloc_props.contiguous()) { auto v = get_view< ST*,HD>(); auto v_src = src.get_view(); - if (src_alloc_size==tgt_alloc_size) { - Kokkos::deep_copy(v,v_src); - } else { - Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { v(idx) = v_src(idx); }); - } } else { auto v = get_strided_view< ST*,HD>(); auto v_src = src.get_strided_view(); - if (src_alloc_size==tgt_alloc_size) { - Kokkos::deep_copy(v,v_src); - } else { - Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { v(idx) = v_src(idx); }); - } } } break; case 2: { - auto v = get_view< ST**,HD>(); - auto v_src = src.get_view(); - if (src_alloc_size==tgt_alloc_size) { - Kokkos::deep_copy(v,v_src); - } else { + if (src_alloc_props.contiguous() and tgt_alloc_props.contiguous()) { + auto v = get_view< ST**,HD>(); + auto v_src = src.get_view(); + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { + int i,j; + unflatten_idx(idx,ext,i,j); + v(i,j) = v_src(i,j); + }); + } + else { + auto v = get_strided_view< ST**,HD>(); + auto v_src = src.get_strided_view(); Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { int i,j; unflatten_idx(idx,ext,i,j); @@ -362,11 +359,17 @@ deep_copy_impl (const Field& src) { break; case 3: { - auto v = get_view< ST***,HD>(); - auto v_src = src.get_view(); - if (src_alloc_size==tgt_alloc_size) { - Kokkos::deep_copy(v,v_src); + if (src_alloc_props.contiguous() and tgt_alloc_props.contiguous()) { + auto v = get_view< ST***,HD>(); + auto v_src = src.get_view(); + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { + int i,j,k; + unflatten_idx(idx,ext,i,j,k); + v(i,j,k) = v_src(i,j,k); + }); } else { + auto v = get_strided_view< ST***,HD>(); + auto v_src = src.get_strided_view(); Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { int i,j,k; unflatten_idx(idx,ext,i,j,k); @@ -377,11 +380,17 @@ deep_copy_impl (const Field& src) { break; case 4: { - auto v = get_view< ST****,HD>(); - auto v_src = src.get_view(); - if (src_alloc_size==tgt_alloc_size) { - Kokkos::deep_copy(v,v_src); + if (src_alloc_props.contiguous() and tgt_alloc_props.contiguous()) { + auto v = get_view< ST****,HD>(); + auto v_src = src.get_view(); + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { + int i,j,k,l; + unflatten_idx(idx,ext,i,j,k,l); + v(i,j,k,l) = v_src(i,j,k,l); + }); } else { + auto v = get_strided_view< ST****,HD>(); + auto v_src = src.get_strided_view(); Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { int i,j,k,l; unflatten_idx(idx,ext,i,j,k,l); @@ -392,11 +401,17 @@ deep_copy_impl (const Field& src) { break; case 5: { - auto v = get_view< ST*****,HD>(); - auto v_src = src.get_view(); - if (src_alloc_size==tgt_alloc_size) { - Kokkos::deep_copy(v,v_src); + if (src_alloc_props.contiguous() and tgt_alloc_props.contiguous()) { + auto v = get_view< ST*****,HD>(); + auto v_src = src.get_view(); + Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { + int i,j,k,l,m; + unflatten_idx(idx,ext,i,j,k,l,m); + v(i,j,k,l,m) = v_src(i,j,k,l,m); + }); } else { + auto v = get_view< ST*****,HD>(); + auto v_src = src.get_view(); Kokkos::parallel_for(policy,KOKKOS_LAMBDA(const int idx) { int i,j,k,l,m; unflatten_idx(idx,ext,i,j,k,l,m); @@ -439,32 +454,57 @@ void Field::deep_copy_impl (const ST value) { break; case 2: { - auto v = get_view(); - Kokkos::deep_copy(v,value); + if (m_header->get_alloc_properties().contiguous()) { + auto v = get_view(); + Kokkos::deep_copy(v,value); + } else { + auto v = get_strided_view(); + Kokkos::deep_copy(v,value); + } } break; case 3: { - auto v = get_view(); - Kokkos::deep_copy(v,value); + if (m_header->get_alloc_properties().contiguous()) { + auto v = get_view(); + Kokkos::deep_copy(v,value); + } else { + auto v = get_strided_view(); + Kokkos::deep_copy(v,value); + } } break; case 4: { - auto v = get_view(); - Kokkos::deep_copy(v,value); + if (m_header->get_alloc_properties().contiguous()) { + auto v = get_view(); + Kokkos::deep_copy(v,value); + } else { + auto v = get_strided_view(); + Kokkos::deep_copy(v,value); + } } break; case 5: { - auto v = get_view(); - Kokkos::deep_copy(v,value); + if (m_header->get_alloc_properties().contiguous()) { + auto v = get_view(); + Kokkos::deep_copy(v,value); + } else { + auto v = get_strided_view(); + Kokkos::deep_copy(v,value); + } } break; case 6: { - auto v = get_view(); - Kokkos::deep_copy(v,value); + if (m_header->get_alloc_properties().contiguous()) { + auto v = get_view(); + Kokkos::deep_copy(v,value); + } else { + auto v = get_strided_view(); + Kokkos::deep_copy(v,value); + } } break; default: diff --git a/components/eamxx/src/share/tests/subfield_tests.cpp b/components/eamxx/src/share/tests/subfield_tests.cpp index f40041a7956..97fe20318c5 100644 --- a/components/eamxx/src/share/tests/subfield_tests.cpp +++ b/components/eamxx/src/share/tests/subfield_tests.cpp @@ -300,6 +300,42 @@ TEST_CASE("field", "") { } } } + SECTION("Subfield deep_copy") { + // ============== + /* Rank-3 view */ + // ============== + std::vector t3 = {COL, CMP, LEV}; + std::vector d3 = {5, 10, 2}; + FieldIdentifier fid3("3d", {t3, d3}, m / s, "some_grid"); + + Field f3a(fid3); + f3a.allocate_view(); + randomize(f3a, engine, pdf); + + Field f3b(fid3); + f3b.allocate_view(); + randomize(f3b, engine, pdf); + + const int idim = 1; + const int sl_beg = 3; + const int sl_end = 6; + + auto sfa = f3a.subfield(idim, sl_beg, sl_end); + auto sv_a = sfa.get_strided_view(); + + auto sfb = f3b.subfield(idim, sl_beg, sl_end); + auto sv_b = sfb.get_strided_view(); + + sfb.deep_copy(sfa); + + for (int i = 0; i < d3[0]; i++) { + for (int j = sl_beg; j < sl_end; j++) { + for (int k = 0; k < d3[2]; k++) { + REQUIRE(sv_a(i, j - sl_beg, k) == sv_b(i, j - sl_beg, k)); + } + } + } + } } // Dynamic Subfields SECTION("dynamic_subfield") {