diff --git a/src/pcms/capi/client.cpp b/src/pcms/capi/client.cpp index 414579bf..60980e31 100644 --- a/src/pcms/capi/client.cpp +++ b/src/pcms/capi/client.cpp @@ -78,7 +78,7 @@ class EmptyFieldFactory : public FieldFactory throw pcms_error( "EmptyFieldFactory::CreateField: requires SimpleFieldData"); } - if (fd->GetDOFAccess()->GetOwnedValues().size() != + if (fd->GetDOFAccess()->GetValues().size() != detail::ExpectedFlatFieldDataSize(*layout_)) { throw pcms_error( "EmptyFieldFactory::CreateField: field data size does not match " diff --git a/src/pcms/coupler/field_communicator.hpp b/src/pcms/coupler/field_communicator.hpp index 8dc8de53..cba19f47 100644 --- a/src/pcms/coupler/field_communicator.hpp +++ b/src/pcms/coupler/field_communicator.hpp @@ -54,6 +54,11 @@ class FieldCommunicator auto data = comm_.Recv(redev::Mode::Synchronous); serializer_.Deserialize(field_, make_const_array_view(data), layout_comm_.GetPermutationArray()); + // Deserialize writes owned DOF holders only, so ghost copies still hold + // their previous values. Anything that evaluates the field reads ghosts + // (an owned element can have vertices another rank owns), so bring them up + // to date. + field_.Synchronize(); } private: diff --git a/src/pcms/coupler/field_exchange_planner.cpp b/src/pcms/coupler/field_exchange_planner.cpp index 88f60d2d..50efee06 100644 --- a/src/pcms/coupler/field_exchange_planner.cpp +++ b/src/pcms/coupler/field_exchange_planner.cpp @@ -29,20 +29,6 @@ struct OutMsg redev::LOs offset; }; -// Returns the mesh entity dimension for DOF local_index based on the -// entity offsets array. ent_offsets[d]..ent_offsets[d+1] is the range of -// DOF indices belonging to mesh entity dimension d. -static int GetMeshEntityDim(LO local_index, const EntOffsetsArray& ent_offsets) -{ - for (int d = 0; d < ent_offsets_len - 1; ++d) { - if (local_index >= static_cast(ent_offsets[d]) && - local_index < static_cast(ent_offsets[d + 1])) { - return d; - } - } - return ent_offsets_len - 1; -} - static size_t GetMessageBlockIndex( LO permutation_entry, Rank1View offsets) { @@ -189,6 +175,7 @@ static ReversePartitionMap2 BuildReversePartitionMap( { PCMS_FUNCTION_TIMER; auto owned = layout.GetOwnedHost(); + auto owned_indices = layout.GetOwnedIndicesHost(); auto class_dims = layout.GetDOFHolderClassificationDimensionsHost(); auto class_ids = layout.GetDOFHolderClassificationIdsHost(); auto coords = layout.GetDOFHolderCoordinates().GetValues(); @@ -209,28 +196,33 @@ static ReversePartitionMap2 BuildReversePartitionMap( Kokkos::create_mirror_view_and_copy(HostMemorySpace(), coords_device); ReversePartitionMap2 reverse_partition; - LO n = static_cast(owned.extent(0)); + const LO n_owned = static_cast(coords.extent(0)); std::array coord{}; auto overlap_mask_view = overlap_mask.GetMask(layout); - for (LO local_index = 0; local_index < n; ++local_index) { + for (LO owned_index = 0; owned_index < n_owned; ++owned_index) { + // Skip holders not owned by this rank. This is required for xgc, whose + // layout marks holders as not owned without having ghosts. The ownership + // mask is indexed with ghosts. + const LO local_index = + owned_indices.size() == 0 ? owned_index : owned_indices(owned_index); if (!owned(local_index)) continue; - if (!overlap_mask_view[local_index]) + if (!overlap_mask_view[owned_index]) continue; for (int d = 0; d < mesh_dim; ++d) - coord[d] = coords_host(local_index, d); + coord[d] = coords_host(owned_index, d); for (int d = mesh_dim; d < 3; ++d) coord[d] = 0.0; - int mesh_ent_dim = GetMeshEntityDim(local_index, ent_offsets); - LO class_dim = class_dims[local_index]; - LO class_id = class_ids[local_index]; + int mesh_ent_dim = GetMeshEntityDim(owned_index, ent_offsets); + LO class_dim = class_dims[owned_index]; + LO class_id = class_ids[owned_index]; auto dr = std::visit(GetRank{class_id, class_dim, coord}, partition); - reverse_partition[dr].indices.emplace_back(local_index); + reverse_partition[dr].indices.emplace_back(owned_index); for (size_t e = static_cast(mesh_ent_dim) + 1; e < ent_offsets_len; ++e) { @@ -296,28 +288,25 @@ void GenericFieldExchangePlanner::FillGidMessage( plan.msg_size + header_size); auto gids = layout.GetGidsHost(); - auto owned = layout.GetOwnedHost(); auto ent_offsets = layout.GetEntOffsets(); auto offsets = Rank1View( plan.offsets.data(), plan.offsets.size()); std::vector per_rank_offsets(plan.dest_ranks.size()); - for (LO local_index = 0; local_index < static_cast(gids.size()); - ++local_index) { - LO perm_index = plan.permutation[local_index]; + for (LO owned_index = 0; owned_index < static_cast(gids.size()); + ++owned_index) { + LO perm_index = plan.permutation[owned_index]; // Owned holders outside the overlap region carry the sentinel and have no // slot in the message. if (perm_index < 0) continue; - // A holder with a valid permutation slot must be owned. - PCMS_ALWAYS_ASSERT(owned[local_index]); auto block_index = GetMessageBlockIndex(perm_index, offsets); const auto gid_index = perm_index + static_cast((block_index + 1) * ent_offsets_len); - gid_message(gid_index) = gids(local_index); + gid_message(gid_index) = gids(owned_index); - int mesh_ent_dim = GetMeshEntityDim(local_index, ent_offsets); + int mesh_ent_dim = GetMeshEntityDim(owned_index, ent_offsets); for (size_t e = static_cast(mesh_ent_dim) + 1; e < ent_offsets_len; ++e) { per_rank_offsets[block_index][e] += 1; diff --git a/src/pcms/coupler/field_layout_communicator.cpp b/src/pcms/coupler/field_layout_communicator.cpp index f20b5394..103fe32e 100644 --- a/src/pcms/coupler/field_layout_communicator.cpp +++ b/src/pcms/coupler/field_layout_communicator.cpp @@ -24,9 +24,9 @@ FieldLayoutCommunicator::FieldLayoutCommunicator( name_(name), redev_(redev), planner_(std::move(planner)), - overlap_mask_(overlap_mask ? std::make_unique(*overlap_mask) - : std::make_unique( - layout.GetGidsHost().size())), + overlap_mask_(overlap_mask + ? std::make_unique(*overlap_mask) + : std::make_unique(layout.GetNumDofHolder())), own_mpi_comm_(own_mpi_comm) { gid_comm_ = channel.CreateComm(name_ + "_gids", mpi_comm_); @@ -111,8 +111,8 @@ void FieldLayoutCommunicator::UpdateLayoutNull() PCMS_FUNCTION_TIMER; // This rank exchanges nothing, but it still unpacks on receive (backends // that keep replicas synchronize them there), so every holder is unmapped. - plan_.permutation = MakeUnmappedPermutation( - static_cast(layout_.GetNumOwnedDofHolder())); + plan_.permutation = + MakeUnmappedPermutation(static_cast(layout_.GetNumDofHolder())); if (redev_.GetProcessType() == redev::ProcessType::Client) { channel_.BeginSendCommunicationPhase(); channel_.EndSendCommunicationPhase(); diff --git a/src/pcms/coupler/field_serializer.h b/src/pcms/coupler/field_serializer.h index 883b6141..a408a425 100644 --- a/src/pcms/coupler/field_serializer.h +++ b/src/pcms/coupler/field_serializer.h @@ -13,7 +13,7 @@ namespace pcms // Packs a field's owned values into, and unpacks them from, an exchange // buffer, through the field's DOFAccess. Any field with DOFAccess is // exchanged this way. Backends that need extra work on receive (e.g., -// synchronizing replicas) do it in DOFAccess::SetOwnedValues. +// synchronizing replicas) do it in DOFAccess::SetValues. template class FieldSerializer { @@ -21,6 +21,7 @@ class FieldSerializer int Serialize(const Field& field, Rank1View buffer, Rank1View permutation) const { + // Only owned (rank-exclusive) DOF holders are serialized. const auto data = field.GetDOFHolderDataHost().GetValues(); if (buffer.size() > 0) { const LO num_dof = static_cast(data.extent(0)); @@ -40,15 +41,14 @@ class FieldSerializer Rank1View permutation) const { const auto& layout = field.GetLayout(); - const LO num_dof = layout.GetNumOwnedDofHolder(); + const LO num_dof = layout.GetNumDofHolder(); const LO num_comp = layout.GetNumComponents(); - // We initialize with local data so that - // data that doesn't get received from deserialized - // doesn't get zeroed out. Long term FIXME + // We initialize with the current data so that data that doesn't get + // received from deserialize doesn't get zeroed out. Long term FIXME // is to make Get/SetData operate over only part of the // arrays const auto current = field.GetDOFHolderDataHost().GetValues(); - Kokkos::View sorted("sorted", layout.OwnedSize()); + Kokkos::View sorted("sorted", layout.Size()); for (LO i = 0; i < num_dof; ++i) { for (LO c = 0; c < num_comp; ++c) { sorted[i * num_comp + c] = current(i, c); diff --git a/src/pcms/coupler/overlap_mask.h b/src/pcms/coupler/overlap_mask.h index 1a917abb..825f6675 100644 --- a/src/pcms/coupler/overlap_mask.h +++ b/src/pcms/coupler/overlap_mask.h @@ -15,6 +15,9 @@ struct OverlapMask Kokkos::View is_overlap_; std::function in_overlap_func_; + // Masks are indexed by owned DOF holder, in owned order, matching the + // owned-indexed exchange plan. + // Default constructor: all DOF holders are in overlap OverlapMask(size_t size) : is_overlap_("overlap_info", size) { @@ -30,7 +33,7 @@ struct OverlapMask } // Construct from a precomputed per-DOF-holder host mask (e.g. an MFEM - // subdomain selected by element attribute). Indexed by local DOF holder. + // subdomain selected by element attribute), one entry per owned DOF holder. OverlapMask(size_t size, Kokkos::View is_overlap_host) : is_overlap_("overlap_info", size) { @@ -38,7 +41,7 @@ struct OverlapMask Kokkos::deep_copy(is_overlap_, is_overlap_host); } - // Construct from Omega_h host array + // Construct from Omega_h host array, one entry per owned DOF holder OverlapMask(size_t size, Omega_h::HostRead is_overlap_host) : is_overlap_("overlap_info", size) { @@ -47,7 +50,7 @@ struct OverlapMask } } - // Construct from Omega_h device array + // Construct from Omega_h device array, one entry per owned DOF holder OverlapMask(size_t size, Omega_h::Read is_overlap_device) : is_overlap_("overlap_info", size) { @@ -57,11 +60,14 @@ struct OverlapMask } } - // Get the mask, evaluating the function against the layout if one was given + // Get the mask, one entry per owned DOF holder, evaluating the function + // against the layout's classification if one was given. Kokkos::View GetMask( const FieldLayout& layout) const { if (!in_overlap_func_) { + PCMS_ALWAYS_ASSERT(is_overlap_.extent(0) == + static_cast(layout.GetNumDofHolder())); return is_overlap_; } auto class_dims = layout.GetDOFHolderClassificationDimensionsHost(); diff --git a/src/pcms/field/CMakeLists.txt b/src/pcms/field/CMakeLists.txt index 496f60e8..00c3f5d7 100644 --- a/src/pcms/field/CMakeLists.txt +++ b/src/pcms/field/CMakeLists.txt @@ -53,6 +53,7 @@ if (PCMS_ENABLE_OMEGA_H) layout/omega_h_entity.h layout/omega_h_lagrange.h evaluator/omega_h_lagrange.h + data/omega_h.h layout/uniform_grid.h uniform_grid_binary_field.h evaluator/uniform_grid.h diff --git a/src/pcms/field/data/mesh_fields.h b/src/pcms/field/data/mesh_fields.h index 0a60b59a..d314611b 100644 --- a/src/pcms/field/data/mesh_fields.h +++ b/src/pcms/field/data/mesh_fields.h @@ -7,6 +7,7 @@ #include "pcms/field/value_view.hpp" #include "pcms/utility/assert.h" #include "pcms/utility/arrays.h" +#include "pcms/utility/omega_h_array_utils.h" #include #include @@ -38,6 +39,32 @@ struct FlattenFunctor } } }; + +// Functor that scatters a flat dof-major buffer back into component-major 2D +// data (LayoutLeft). Inverse of FlattenFunctor; indexes the full destination +// view directly (rather than via Kokkos::subview) because a partial row-range +// subview of a LayoutLeft view is not representable as LayoutLeft/LayoutRight +// (Kokkos falls back to LayoutStride), which Rank2View does not detect. +template +struct UnflattenFunctor +{ + Kokkos::View flat_; + Rank2View data_; + LO row_off_; + size_t nc_; + UnflattenFunctor(Kokkos::View f, + Rank2View d, LO ro, size_t nc) + : flat_(f), data_(d), row_off_(ro), nc_(nc) + { + } + KOKKOS_INLINE_FUNCTION void operator()(LO local) const + { + LO global_dof = row_off_ + local; + for (size_t c = 0; c < nc_; ++c) { + data_(global_dof, c) = flat_(local * nc_ + c); + } + } +}; } // namespace template @@ -46,14 +73,17 @@ class MeshFieldsFieldData , public DOFAccess { public: + // Owned values are the stored data. The MeshFields backend holds the values + // of every DOF holder, owned and ghost, which its evaluator reads: it + // shares the owned values and gets the ghost values from Synchronize. MeshFieldsFieldData(std::shared_ptr layout, ValueBasis basis) : layout_(std::move(layout)), basis_(std::move(basis)), mesh_field_(MakeMeshFieldBackend(*layout_)), - device_data_("meshfields_field_data_device", - static_cast(layout_->GetNumOwnedDofHolder()), - static_cast(layout_->GetNumComponents())) + owned_data_("meshfields_field_data_device", + static_cast(layout_->GetNumDofHolder()), + static_cast(layout_->GetNumComponents())) { if (!mesh_field_) { throw pcms_error( @@ -70,15 +100,48 @@ class MeshFieldsFieldData const DOFAccess* GetDOFAccess() const noexcept override { return this; } DOFAccess* GetDOFAccess() noexcept override { return this; } - Rank2View GetOwnedValues() const override + Rank2View GetValues() const override { - return MakeConstRank2View(device_data_); + return MakeConstRank2View(owned_data_); } - void SetOwnedValues(Rank2View values) override + void SetValues(Rank2View values) override { - CopyDeviceRank2ViewToDeviceView(device_data_, values); - SyncBackend(GetOwnedValues()); + CopyDeviceRank2ViewToDeviceView(owned_data_, values); + if (!layout_->HasGhosts()) { + SyncBackend(MakeConstRank2View(owned_data_)); + return; + } + // Write the owned rows into the backend's values with ghosts, keeping its + // ghost values. + auto local = ReadBackend(); + ScatterRows(MakeConstRank2View(owned_data_), layout_->GetOwnedIndices(), + MakeRank2View(local)); + SyncBackend(MakeConstRank2View(local)); + } + + void Synchronize() override + { + const int nc = layout_->GetNumComponents(); + auto& mesh = layout_->GetMesh(); + const auto nodes_per_dim = layout_->GetNodesPerDim(); + for (int dim = 0; dim <= mesh.dim(); ++dim) { + if (!nodes_per_dim[dim]) { + continue; + } + // The backend's per-entity data is already dof-major, the order + // Omega_h's sync_array works in. + const int width = nodes_per_dim[dim] * nc; + const size_t len = + static_cast(mesh.nents(dim)) * static_cast(width); + Kokkos::View flat("ghost_sync_flat", len); + mesh_field_->GetData(Rank1View(flat.data(), len), + nodes_per_dim[dim], nc, dim); + SynchronizeOmegaHBlock(mesh, dim, width, flat); + mesh_field_->SetData( + Rank1View(flat.data(), len), + nodes_per_dim[dim], nc, dim); + } } std::shared_ptr> GetMeshFieldBackend() const @@ -87,6 +150,37 @@ class MeshFieldsFieldData } private: + // The backend's values of every DOF holder with ghosts, in mesh order. + Kokkos::View ReadBackend() const + { + const auto nodes_per_dim = layout_->GetNodesPerDim(); + const auto nc = static_cast(layout_->GetNumComponents()); + auto& mesh = layout_->GetMesh(); + Kokkos::View local( + "meshfields_local_values", + static_cast(layout_->GetNumDofHolderWithGhosts()), nc); + size_t row_offset = 0; + for (int i = 0; i <= mesh.dim(); ++i) { + if (!nodes_per_dim[i]) { + continue; + } + const size_t num_rows = static_cast(mesh.nents(i)) * + static_cast(nodes_per_dim[i]); + const size_t len = num_rows * nc; + Kokkos::View flat("read_flat", len); + mesh_field_->GetData(Rank1View(flat.data(), len), + nodes_per_dim[i], nc, i); + Kokkos::parallel_for( + "ReadBackendReorder", + Kokkos::RangePolicy( + 0, static_cast(num_rows)), + UnflattenFunctor(flat, MakeRank2View(local), + static_cast(row_offset), nc)); + row_offset += num_rows; + } + return local; + } + // Serialization boundary: meshfields' SetData consumes flat node-major // spans, so the shaped data is explicitly repacked into a flat node-major // staging buffer here — the one place this backend handles flat memory. @@ -95,7 +189,7 @@ class MeshFieldsFieldData auto nodes_per_dim = layout_->GetNodesPerDim(); auto num_components = layout_->GetNumComponents(); auto& mesh = layout_->GetMesh(); - // device_data_ is rank-2 LayoutLeft (component-major), but SetData + // The shaped data may be LayoutLeft (component-major), but SetData // expects a flat dof-major span. size_t row_offset = 0; for (int i = 0; i <= mesh.dim(); ++i) { @@ -120,7 +214,7 @@ class MeshFieldsFieldData std::shared_ptr layout_; ValueBasis basis_; std::shared_ptr> mesh_field_; - Kokkos::View device_data_; + Kokkos::View owned_data_; }; } // namespace pcms diff --git a/src/pcms/field/data/mfem.cpp b/src/pcms/field/data/mfem.cpp index 00405bad..cc5539c9 100644 --- a/src/pcms/field/data/mfem.cpp +++ b/src/pcms/field/data/mfem.cpp @@ -78,15 +78,14 @@ void MFEMVertexFieldData::ScatterHost( gf_.SetFromTrueDofs(true_values); } -Rank2View MFEMVertexFieldData::GetOwnedValues() - const +Rank2View MFEMVertexFieldData::GetValues() const { GatherHost(); DeepCopyMismatchLayouts(device_data_, host_data_); return MakeConstRank2View(device_data_); } -void MFEMVertexFieldData::SetOwnedValues( +void MFEMVertexFieldData::SetValues( Rank2View data) { CopyDeviceRank2ViewToDeviceView(device_data_, data); diff --git a/src/pcms/field/data/mfem.hpp b/src/pcms/field/data/mfem.hpp index be7837a1..b0a4ce34 100644 --- a/src/pcms/field/data/mfem.hpp +++ b/src/pcms/field/data/mfem.hpp @@ -25,8 +25,8 @@ class MFEMVertexFieldData const DOFAccess* GetDOFAccess() const noexcept override { return this; } DOFAccess* GetDOFAccess() noexcept override { return this; } - Rank2View GetOwnedValues() const override; - void SetOwnedValues(Rank2View data) override; + Rank2View GetValues() const override; + void SetValues(Rank2View data) override; private: Rank2View GatherHost() const; diff --git a/src/pcms/field/data/omega_h.h b/src/pcms/field/data/omega_h.h new file mode 100644 index 00000000..32216b8a --- /dev/null +++ b/src/pcms/field/data/omega_h.h @@ -0,0 +1,157 @@ +#ifndef PCMS_OMEGA_H_FIELD_DATA_H +#define PCMS_OMEGA_H_FIELD_DATA_H + +#include "../field_data.h" +#include "../field_layout.h" +#include "pcms/discretization/discretization/omega_h.hpp" +#include "pcms/utility/arrays.h" +#include "pcms/utility/memory_spaces.h" +#include "pcms/utility/omega_h_array_utils.h" +#include +#include +#include +#include +#include +#include + +namespace pcms +{ + +namespace detail +{ +// A tag array (one value per component per entity, entity-major) holding the +// owned values at their entities and zero at the ghosts. Local to this rank. +template +Omega_h::Read OwnedValuesAsTag( + Rank2View owned, + Kokkos::View owned_indices, LO num_entities, + int num_components) +{ + Omega_h::Write tag(num_entities * num_components, TagValue(0)); + Kokkos::parallel_for( + "pcms::OwnedValuesAsTag", + Kokkos::RangePolicy( + 0, static_cast(owned.extent(0))), + KOKKOS_LAMBDA(LO o) { + const LO i = owned_indices(o); + for (int c = 0; c < num_components; ++c) { + tag[i * num_components + c] = static_cast(owned(o, c)); + } + }); + return tag; +} + +// A name for a tag an OmegaHFieldData creates, unique within the process. +inline std::string NextOmegaHFieldDataTagName() +{ + static std::atomic counter{0}; + return "pcms_field_" + std::to_string(counter++); +} +} // namespace detail + +// OmegaHFieldData is the FieldData of an Omega_h layout. It stores the +// owned values in a [dof_holder][component] Kokkos View and also keeps the +// values of every DOF holder, owned and ghost, in an Omega_h tag, for code that +// reads ghosts (e.g. the Omega_h Lagrange evaluator). +// +// It creates a tag of its own on the mesh, named by +// detail::NextOmegaHFieldDataTagName(), and removes it on destruction, so the +// mesh must outlive it. +// SetValues replaces the tag with the owned values (ghosts zero, local to this +// rank); Synchronize (collective) fills the ghosts from their owners. Ghost +// values are correct only after Synchronize. +template +class OmegaHFieldData + : public FieldData + , public DOFAccess +{ +public: + // Value type of the Omega_h tag. + using tag_value_type = OmegaHValue_t; + + OmegaHFieldData(std::shared_ptr layout, ValueBasis basis) + : layout_(std::move(layout)), + basis_(std::move(basis)), + device_data_("omega_h_field_data_device", + static_cast(layout_->GetNumDofHolder()), + static_cast(layout_->GetNumComponents())) + { + const auto* oh = dynamic_cast( + layout_->GetDiscretization().get()); + if (oh == nullptr) { + throw pcms_error("OmegaHFieldData: layout is not on an Omega_h mesh"); + } + tag_dim_ = layout_->GetDOFHolderEntityDim(); + if (tag_dim_ < 0) { + throw pcms_error("OmegaHFieldData: layout DOF holders are not one per " + "mesh entity of a single dimension"); + } + mesh_ = &oh->GetMesh(); + tag_name_ = detail::NextOmegaHFieldDataTagName(); + const int nc = layout_->GetNumComponents(); + mesh_->add_tag( + tag_dim_, tag_name_, nc, + Omega_h::Read(mesh_->nents(tag_dim_) * nc, + tag_value_type(0))); + } + + // The tag holds the name of this object, so it is neither copied nor moved. + OmegaHFieldData(const OmegaHFieldData&) = delete; + OmegaHFieldData& operator=(const OmegaHFieldData&) = delete; + OmegaHFieldData(OmegaHFieldData&&) = delete; + OmegaHFieldData& operator=(OmegaHFieldData&&) = delete; + + ~OmegaHFieldData() noexcept override + { + if (mesh_->has_tag(tag_dim_, tag_name_)) { + mesh_->remove_tag(tag_dim_, tag_name_); + } + } + + FieldValueType GetValueType() const override + { + return ValueTypeOfRank(basis_.Rank()); + } + const ValueBasis& GetValueBasis() const override { return basis_; } + + const DOFAccess* GetDOFAccess() const noexcept override { return this; } + DOFAccess* GetDOFAccess() noexcept override { return this; } + + Rank2View GetValues() const override + { + return MakeConstRank2View(device_data_); + } + + void SetValues(Rank2View values) override + { + CopyDeviceRank2ViewToDeviceView(device_data_, values); + mesh_->set_tag(tag_dim_, tag_name_, + detail::OwnedValuesAsTag( + GetValues(), layout_->GetOwnedIndices(), + mesh_->nents(tag_dim_), layout_->GetNumComponents())); + } + + // Fills the tag's ghost entries with their owners' values. Collective. + void Synchronize() override { mesh_->sync_tag(tag_dim_, tag_name_); } + + // Name of the Omega_h tag holding the values of every DOF holder, owned and + // ghost, on entity dimension layout.GetDOFHolderEntityDim() (one value per + // component per entity). Its owned entries + // are always current; its ghost entries are correct only after Synchronize. + [[nodiscard]] const std::string& GetTagName() const noexcept + { + return tag_name_; + } + +private: + std::shared_ptr layout_; + ValueBasis basis_; + Kokkos::View device_data_; + Omega_h::Mesh* mesh_ = nullptr; + int tag_dim_ = -1; + std::string tag_name_; +}; + +} // namespace pcms + +#endif // PCMS_OMEGA_H_FIELD_DATA_H diff --git a/src/pcms/field/data/point_cloud.cpp b/src/pcms/field/data/point_cloud.cpp index d50e512d..79a36af4 100644 --- a/src/pcms/field/data/point_cloud.cpp +++ b/src/pcms/field/data/point_cloud.cpp @@ -25,12 +25,12 @@ const ValueBasis& PointCloud::GetValueBasis() const return basis_; } -Rank2View PointCloud::GetOwnedValues() const +Rank2View PointCloud::GetValues() const { return MakeConstRank2View(device_data_); } -void PointCloud::SetOwnedValues(Rank2View data) +void PointCloud::SetValues(Rank2View data) { PCMS_FUNCTION_TIMER; CopyDeviceRank2ViewToDeviceView(device_data_, data); diff --git a/src/pcms/field/data/point_cloud.h b/src/pcms/field/data/point_cloud.h index a7006b9e..8eda5f39 100644 --- a/src/pcms/field/data/point_cloud.h +++ b/src/pcms/field/data/point_cloud.h @@ -22,8 +22,8 @@ class PointCloud const DOFAccess* GetDOFAccess() const noexcept override { return this; } DOFAccess* GetDOFAccess() noexcept override { return this; } - Rank2View GetOwnedValues() const override; - void SetOwnedValues(Rank2View data) override; + Rank2View GetValues() const override; + void SetValues(Rank2View data) override; private: std::shared_ptr layout_; diff --git a/src/pcms/field/data/simple.h b/src/pcms/field/data/simple.h index facfd032..dfd2ff26 100644 --- a/src/pcms/field/data/simple.h +++ b/src/pcms/field/data/simple.h @@ -28,7 +28,7 @@ class SimpleFieldData : layout_(std::move(layout)), basis_(std::move(basis)), device_data_("simple_field_data_device", - static_cast(layout_->GetNumOwnedDofHolder()), + static_cast(layout_->GetNumDofHolder()), static_cast(layout_->GetNumComponents())) { } @@ -42,12 +42,12 @@ class SimpleFieldData const DOFAccess* GetDOFAccess() const noexcept override { return this; } DOFAccess* GetDOFAccess() noexcept override { return this; } - Rank2View GetOwnedValues() const override + Rank2View GetValues() const override { return MakeConstRank2View(device_data_); } - void SetOwnedValues(Rank2View values) override + void SetValues(Rank2View values) override { CopyDeviceRank2ViewToDeviceView(device_data_, values); } diff --git a/src/pcms/field/data/xgc.h b/src/pcms/field/data/xgc.h index 9caf89d7..4bb2a09c 100644 --- a/src/pcms/field/data/xgc.h +++ b/src/pcms/field/data/xgc.h @@ -23,7 +23,7 @@ class XGCFieldData // // XGC keeps a full copy of a plane's data on every rank of the plane, and // only the plane's root rank (rank 0 of plane_comm) exchanges it. When - // plane_comm is given, SetOwnedValues broadcasts the root's values to the + // plane_comm is given, SetValues broadcasts the root's values to the // other ranks of the plane, so every rank of plane_comm must call it. XGCFieldData(std::shared_ptr layout, ValueBasis basis, Rank1View data, @@ -59,7 +59,7 @@ class XGCFieldData const DOFAccess* GetDOFAccess() const noexcept override { return this; } DOFAccess* GetDOFAccess() noexcept override { return this; } - Rank2View GetOwnedValues() const override + Rank2View GetValues() const override { EnsureDeviceStaging(); const auto nc = layout_->GetNumComponents(); @@ -70,7 +70,7 @@ class XGCFieldData return MakeConstRank2View(device_data_); } - void SetOwnedValues(Rank2View values) override + void SetValues(Rank2View values) override { EnsureDeviceStaging(); CopyDeviceRank2ViewToDeviceView(device_data_, values); diff --git a/src/pcms/field/evaluator/mesh_fields_backend.h b/src/pcms/field/evaluator/mesh_fields_backend.h index e59025e8..c6478eba 100644 --- a/src/pcms/field/evaluator/mesh_fields_backend.h +++ b/src/pcms/field/evaluator/mesh_fields_backend.h @@ -14,7 +14,8 @@ #include "pcms/utility/assert.h" #include "pcms/utility/arrays.h" #include "pcms/localization/point_search.h" -#include "pcms/field/field.h" // OutOfBoundsMode +#include "pcms/field/field.h" +#include "pcms/field/out_of_bounds_policy.h" namespace pcms { diff --git a/src/pcms/field/evaluator/omega_h_lagrange.h b/src/pcms/field/evaluator/omega_h_lagrange.h index 72bc95db..ce2d47fd 100644 --- a/src/pcms/field/evaluator/omega_h_lagrange.h +++ b/src/pcms/field/evaluator/omega_h_lagrange.h @@ -6,6 +6,7 @@ #include "pcms/field/out_of_bounds_policy.h" #include "pcms/field/point_evaluator.h" #include "pcms/field/field_data.h" +#include "pcms/field/data/omega_h.h" #include "pcms/localization/point_search.h" #include "pcms/utility/assert.h" #include "pcms/utility/arrays.h" @@ -199,12 +200,33 @@ inline std::variant MakeSearch( // PointEvaluator // --------------------------------------------------------------------------- +namespace detail +{ +// DOF-holder values read from an Omega_h tag: one value per component per +// mesh entity, entity-major. +template +struct OmegaHLagrangeTagValues +{ + Omega_h::Read values; + int num_components; + KOKKOS_INLINE_FUNCTION T operator()(LO holder, int c) const + { + return static_cast(values[holder * num_components + c]); + } +}; +} // namespace detail + // OmegaHLagrangePointEvaluator implements PointEvaluator for simplex // meshes backed by Omega_h. Localization results (element IDs, barycentric // coordinates) are computed once at construction and cached for repeated // Evaluate calls. // // Output shape: [num_query_points][num_components]. +// +// The values are read from the field's Omega_h tag +// (OmegaHFieldData::GetTagName). Its owned entries are always current; its +// ghost entries, which elements on the partition boundary read, are correct +// only after the field's last write was followed by Synchronize. template > @@ -235,7 +257,14 @@ class OmegaHLagrangePointEvaluator : public PointEvaluator { PCMS_FUNCTION_TIMER; const auto values = out.GetValues(); - const auto dof_data = field.GetDOFHolderData().GetValues(); + // Elements read the values of all their DOF holders, ghosts included, so + // read the field's Omega_h tag rather than the owned values. + const auto& data = dynamic_cast&>(field.GetData()); + using TagValue = typename OmegaHFieldData::tag_value_type; + const auto dof_data = detail::OmegaHLagrangeTagValues{ + layout_->GetMesh().get_array(layout_->GetDOFHolderEntityDim(), + data.GetTagName()), + layout_->GetNumComponents()}; LO n_valid = static_cast(hint_.elem_ids.size()); int n_comp = layout_->GetNumComponents(); diff --git a/src/pcms/field/evaluator/uniform_grid.h b/src/pcms/field/evaluator/uniform_grid.h index 35fdd293..9c364ba3 100644 --- a/src/pcms/field/evaluator/uniform_grid.h +++ b/src/pcms/field/evaluator/uniform_grid.h @@ -131,7 +131,7 @@ class UniformGridPointEvaluator : public PointEvaluator int n_comp = layout_->GetNumComponents(); PCMS_ALWAYS_ASSERT(dof_data.size() == - static_cast(layout_->GetNumOwnedDofHolder() * + static_cast(layout_->GetNumDofHolder() * layout_->GetNumComponents())); auto cell_indices = hint_.cell_indices_; @@ -208,13 +208,13 @@ class UniformGridPointEvaluator : public PointEvaluator cell_indices_interp(i, d) = cell_dim_indices(i, d); }); - RealMatView values_interp( - "values_interp", static_cast(layout_->GetNumOwnedDofHolder()), - static_cast(n_comp)); + RealMatView values_interp("values_interp", + static_cast(layout_->GetNumDofHolder()), + static_cast(n_comp)); Kokkos::parallel_for( "copy_dof_data_to_values_interp", Kokkos::MDRangePolicy>( - {0, 0}, {static_cast(layout_->GetNumOwnedDofHolder()), n_comp}), + {0, 0}, {static_cast(layout_->GetNumDofHolder()), n_comp}), KOKKOS_LAMBDA(const LO i, const int c) { values_interp(i, c) = dof_data(i, c); }); diff --git a/src/pcms/field/field.h b/src/pcms/field/field.h index ab70297b..e8423143 100644 --- a/src/pcms/field/field.h +++ b/src/pcms/field/field.h @@ -43,52 +43,59 @@ class Field [[nodiscard]] const std::string& GetName() const noexcept { return name_; } // The DOF-holder accessors below read and write the values of the DOF - // holders this rank owns through the data's DOFAccess, and throw if the data - // has none (e.g. a field of an evaluation-only source). Writes synchronize - // the backend's other copies and are collective; see DOFAccess. + // holders this rank owns, shaped [GetNumDofHolder()][component] in owned + // order, through the data's DOFAccess, and throw if the data has none (e.g. + // a field of an evaluation-only source). Writes synchronize the backend's + // other copies and are collective; see DOFAccess. Ghost copies are brought + // up to date by Synchronize. // Copies to the host unless device memory is host memory. The returned view // is valid until the next call or until the field is written or destroyed. [[nodiscard]] ValueView GetDOFHolderDataHost() const { return ValueView(data_->GetValueBasis(), - OwnedValuesHost()); + ValuesHost()); } [[nodiscard]] ValueView GetDOFHolderData() const { return ValueView( - data_->GetValueBasis(), RequireDOFAccess().GetOwnedValues()); + data_->GetValueBasis(), RequireDOFAccess().GetValues()); } void SetDOFHolderDataHost(ValueView v) { CheckWriteTag(v.GetBasis()); CheckWriteExtents(v.GetValues().extent(0), v.GetValues().extent(1)); - SetOwnedValuesFromHost(v.GetValues()); + SetValuesFromHost(v.GetValues()); } void SetDOFHolderData(ValueView v) { CheckWriteTag(v.GetBasis()); CheckWriteExtents(v.GetValues().extent(0), v.GetValues().extent(1)); - RequireDOFAccess().SetOwnedValues(v.GetValues()); + RequireDOFAccess().SetValues(v.GetValues()); } // unchecked is needed for optimization in some places, but should generally // not be used unless you know what you are doing void SetDOFHolderDataUncheckedHost(Rank2View v) { - SetOwnedValuesFromHost(v); + SetValuesFromHost(v); } // unchecked is needed for optimization in some places, but should generally // not be used unless you know what you are doing void SetDOFHolderDataUnchecked(Rank2View v) { - RequireDOFAccess().SetOwnedValues(v); + RequireDOFAccess().SetValues(v); } + // Bring the backend's copies of DOF holders other ranks own (e.g. Omega_h + // ghosts) up to date with their owners' values. Collective; a no-op for + // non-distributed fields. + void Synchronize() { data_->Synchronize(); } + protected: Field(std::string name, std::shared_ptr layout, std::unique_ptr> data) @@ -108,13 +115,13 @@ class Field void CheckWriteExtents(size_t num_holders, size_t num_components) const { - if (num_holders != static_cast(layout_->GetNumOwnedDofHolder()) || + if (num_holders != static_cast(layout_->GetNumDofHolder()) || num_components != static_cast(layout_->GetNumComponents())) { throw std::invalid_argument( "Field::SetDOFHolderData: the written view is shaped [" + std::to_string(num_holders) + "][" + std::to_string(num_components) + "] but this field's layout is [" + - std::to_string(layout_->GetNumOwnedDofHolder()) + "][" + + std::to_string(layout_->GetNumDofHolder()) + "][" + std::to_string(layout_->GetNumComponents()) + "]"); } } @@ -135,9 +142,9 @@ class Field return const_cast&>(std::as_const(*this).RequireDOFAccess()); } - Rank2View OwnedValuesHost() const + Rank2View ValuesHost() const { - const auto values = RequireDOFAccess().GetOwnedValues(); + const auto values = RequireDOFAccess().GetValues(); if constexpr (std::is_same_v) { return values; } else { @@ -155,16 +162,16 @@ class Field } } - void SetOwnedValuesFromHost(Rank2View values) + void SetValuesFromHost(Rank2View values) { if constexpr (std::is_same_v) { - RequireDOFAccess().SetOwnedValues(values); + RequireDOFAccess().SetValues(values); } else { Kokkos::View staged( Kokkos::view_alloc(Kokkos::WithoutInitializing, "owned_values_staged"), values.extent(0), values.extent(1)); CopyHostRank2ViewToDeviceView(staged, values); - RequireDOFAccess().SetOwnedValues(MakeConstRank2View(staged)); + RequireDOFAccess().SetValues(MakeConstRank2View(staged)); } } diff --git a/src/pcms/field/field_data.h b/src/pcms/field/field_data.h index bed91f93..7d1efc18 100644 --- a/src/pcms/field/field_data.h +++ b/src/pcms/field/field_data.h @@ -19,12 +19,13 @@ class DOFAccess { public: // Valid until the field data is next written or destroyed. - virtual Rank2View GetOwnedValues() const = 0; + virtual Rank2View GetValues() const = 0; - // Writes the owned values, then brings any other copies the backend keeps - // (ghosts, replicas on other ranks) up to date. Collective over the ranks - // that share those copies. - virtual void SetOwnedValues(Rank2View values) = 0; + // Writes the owned values, then brings any internal copies the backend + // keeps (e.g. replicas on other ranks) up to date. Ghost copies keep their + // values until FieldData::Synchronize. Collective over the ranks that + // share those copies. + virtual void SetValues(Rank2View values) = 0; virtual ~DOFAccess() noexcept = default; }; @@ -43,6 +44,11 @@ class FieldData virtual const DOFAccess* GetDOFAccess() const noexcept { return nullptr; } virtual DOFAccess* GetDOFAccess() noexcept { return nullptr; } + // Bring any copies this rank keeps of DOF holders other ranks own (e.g. the + // ghost copies of an Omega_h mesh) up to date with their owners' values. + // Collective; a no-op for backends that keep no such copies. + virtual void Synchronize() {} + virtual ~FieldData() noexcept = default; }; diff --git a/src/pcms/field/field_factory.h b/src/pcms/field/field_factory.h index 93667263..efae413b 100644 --- a/src/pcms/field/field_factory.h +++ b/src/pcms/field/field_factory.h @@ -18,7 +18,7 @@ namespace detail inline size_t ExpectedFlatFieldDataSize(const FieldLayout& layout) { - return static_cast(layout.GetNumOwnedDofHolder()) * + return static_cast(layout.GetNumDofHolder()) * static_cast(layout.GetNumComponents()); } diff --git a/src/pcms/field/field_layout.cpp b/src/pcms/field/field_layout.cpp index 6e785f9e..1463227b 100644 --- a/src/pcms/field/field_layout.cpp +++ b/src/pcms/field/field_layout.cpp @@ -5,10 +5,16 @@ namespace pcms { -void FieldLayout::BuildGlobalToLocalPermutation() +void FieldLayout::BuildGlobalToLocalPermutation( + Kokkos::View local_gids, + const EntOffsetsArray& local_ent_offsets) { - const auto permutation = - BuildSortedGidPermutation(GetGidsHost(), GetEntOffsets()); + // The sort runs on the host; the GIDs are copied there only for it. + const auto gids_host = + Kokkos::create_mirror_view_and_copy(HostMemorySpace(), local_gids); + const auto permutation = BuildSortedGidPermutation( + GlobalIDView(gids_host.data(), gids_host.size()), + local_ent_offsets); global_to_local_host_ = Kokkos::View( "global_to_local_host", permutation.size()); for (size_t i = 0; i < permutation.size(); ++i) diff --git a/src/pcms/field/field_layout.h b/src/pcms/field/field_layout.h index 1c847734..88a0ca0e 100644 --- a/src/pcms/field/field_layout.h +++ b/src/pcms/field/field_layout.h @@ -18,6 +18,45 @@ using EntOffsetsArray = std::array; using ReversePartitionMap = std::map>; +// Returns the mesh entity dimension for the DOF holder at local_index, based on +// the entity-offset array: ent_offsets[d]..ent_offsets[d+1] is the range of DOF +// indices belonging to mesh entity dimension d. +inline int GetMeshEntityDim(LO local_index, const EntOffsetsArray& ent_offsets) +{ + for (int d = 0; d < ent_offsets_len - 1; ++d) { + if (local_index >= static_cast(ent_offsets[d]) && + local_index < static_cast(ent_offsets[d + 1])) { + return d; + } + } + return ent_offsets_len - 1; +} + +// FieldLayout describes how a field's DOF holders are laid out on this rank. +// +// The accessors describe the DOF holders this rank owns (rank-exclusive), +// compact and in owned order. A distributed layout also holds ghost copies of +// DOF holders other ranks own, for internal use (field storage, ghost +// synchronization, element assembly): together they are the DOF holders with +// ghosts, in mesh order. Only their ownership and the owned indices are exposed +// here. For a layout without ghosts the two coincide (!HasGhosts()). +// Entity offsets over the owned DOF holders, in owned order, from the entity +// offsets of the DOF holders with ghosts and the owned indices. +inline EntOffsetsArray OwnedEntOffsets( + const EntOffsetsArray& local_offsets, + Kokkos::View owned_indices) +{ + EntOffsetsArray offsets{}; + offsets.fill(0); + for (size_t o = 0; o < owned_indices.size(); ++o) { + const int d = GetMeshEntityDim(owned_indices(o), local_offsets); + for (size_t e = static_cast(d) + 1; e < ent_offsets_len; ++e) { + offsets[e] += 1; + } + } + return offsets; +} + class FieldLayout { public: @@ -33,19 +72,47 @@ class FieldLayout // number of components int virtual GetNumComponents() const = 0; - // nodes for standard lagrange FEM - LO virtual GetNumOwnedDofHolder() const = 0; + // number of DOF holders this rank owns + LO virtual GetNumDofHolder() const = 0; + GO virtual GetNumGlobalDofHolder() const = 0; - // size of buffer that needs to be allocated to represent the field - // # components * NumDOFHolder - LO OwnedSize() const { return GetNumComponents() * GetNumOwnedDofHolder(); }; + // size of buffer needed to hold the owned coefficients + // # components * NumDofHolder + LO Size() const { return GetNumComponents() * GetNumDofHolder(); } + GO GlobalSize() const { return GetNumComponents() * GetNumGlobalDofHolder(); }; + // True when the layout holds ghost DOF holders, i.e. owned order differs + // from the order with ghosts. Layouts that hold ghosts override it. + [[nodiscard]] virtual bool HasGhosts() const { return false; } + + // returns true if the field layout is distributed: on several ranks it may + // hold ghost DOF holders in addition to the owned ones (see HasGhosts) + [[nodiscard]] virtual bool IsDistributed() const = 0; + + // Per DOF holder, ghosts included: whether this rank owns it. virtual Rank1View GetOwnedHost() const = 0; + + // Maps owned index (0..GetNumDofHolder()-1) to its index among the DOF + // holders with ghosts. An empty view means the identity map (no ghosts). A + // distributed layout on a rank that owns nothing also returns an empty view, + // so callers that need to know whether the layout has no ghosts must ask + // `!HasGhosts()`, not test for empty. + virtual Kokkos::View GetOwnedIndicesHost() const + { + return {}; + } + + // Device analogue of GetOwnedIndicesHost(). + virtual Kokkos::View GetOwnedIndices() const + { + return {}; + } + virtual GlobalIDView GetGidsHost() const = 0; // Maps each local DOF holder to its contiguous active index, ordered by GID @@ -58,10 +125,6 @@ class FieldLayout Kokkos::View GetGlobalToLocalPermutation() const; - // returns true if the field layout is distributed - // if the field layout is distributed, the owned and global dofs are the same - [[nodiscard]] virtual bool IsDistributed() const = 0; - virtual EntOffsetsArray GetEntOffsets() const = 0; [[nodiscard]] const std::shared_ptr& @@ -78,6 +141,13 @@ class FieldLayout virtual int GetDimension() const = 0; + // Entity dimension of the DOF holders (0 = vertex, ..., mesh dim = element) + // when the layout has exactly one DOF holder per local mesh entity of that + // dimension, ghosts included; used for ghost synchronization and to pick + // mesh-based localization. -1 otherwise (e.g. non-mesh or multi-dimension + // layouts). + virtual int GetDOFHolderEntityDim() const { return -1; } + virtual Rank1View GetDOFHolderClassificationDimensionsHost() const = 0; @@ -87,12 +157,15 @@ class FieldLayout virtual ~FieldLayout() noexcept = default; protected: - // Builds the global-to-local permutation from GetGidsHost()/GetEntOffsets(). + // Builds the global-to-local permutation from the GIDs and entity offsets of + // the DOF holders with ghosts (the permutation is indexed with ghosts). // Derived layouts that support GetGlobalToLocalPermutation() must call this // from their constructor, once their GIDs and entity offsets are available; // building eagerly keeps first access free of the data race that lazy // construction of the mutable cache would introduce. - void BuildGlobalToLocalPermutation(); + void BuildGlobalToLocalPermutation( + Kokkos::View local_gids, + const EntOffsetsArray& local_ent_offsets); void SetCoordinateSystem( std::shared_ptr coordinate_system) @@ -111,5 +184,101 @@ class FieldLayout Kokkos::View global_to_local_; }; +// The geometric classification (dimension and id) of the owned DOF holders, +// as host arrays in owned order. Converted and gathered on the device from the +// per-holder (ghosts included) device arrays; only the owned entries are +// copied to the host. +template +void BuildOwnedClassificationHost( + const DimsView& class_dims, const IdsView& class_ids, + Kokkos::View owned_indices, + Kokkos::View& owned_dims_host, + Kokkos::View& owned_ids_host) +{ + using exec = DeviceMemorySpace::execution_space; + const auto num_owned = owned_indices.extent(0); + Kokkos::View owned_dims("owned_class_dims", + num_owned); + Kokkos::View owned_ids("owned_class_ids", num_owned); + Kokkos::parallel_for( + "pcms::GatherOwnedClassification", Kokkos::RangePolicy(0, num_owned), + KOKKOS_LAMBDA(LO o) { + const LO i = owned_indices(o); + owned_dims(o) = + static_cast(static_cast(class_dims[i])); + owned_ids(o) = static_cast(class_ids[i]); + }); + owned_dims_host = + Kokkos::create_mirror_view_and_copy(HostMemorySpace(), owned_dims); + owned_ids_host = + Kokkos::create_mirror_view_and_copy(HostMemorySpace(), owned_ids); +} + +// Compact "owned" (rank-exclusive) views derived from a layout's arrays with +// ghosts. +struct OwnedLayoutData +{ + LO num_owned = 0; + Kokkos::View owned_indices_host; + Kokkos::View owned_gids_host; + Kokkos::View owned_coords_2d; + Kokkos::View owned_indices; +}; + +// Built on the device from the device-resident ownership mask, GIDs and +// coordinates; only the owned indices and owned GIDs are copied to the host, +// for their host-side users. +template +OwnedLayoutData BuildOwnedLayoutData( + Kokkos::View owned, + Kokkos::View gids, const CoordsView& coords, + int dim) +{ + using exec = DeviceMemorySpace::execution_space; + const LO n = static_cast(owned.extent(0)); + OwnedLayoutData out; + + Kokkos::parallel_reduce( + "pcms::CountOwned", Kokkos::RangePolicy(0, n), + KOKKOS_LAMBDA(LO i, LO & count) { count += owned(i) ? 1 : 0; }, + out.num_owned); + + // Owned holders in increasing order: an exclusive scan over the mask gives + // each owned holder its owned index. + out.owned_indices = + Kokkos::View("owned_indices", out.num_owned); + const auto owned_indices = out.owned_indices; + Kokkos::parallel_scan( + "pcms::OwnedIndices", Kokkos::RangePolicy(0, n), + KOKKOS_LAMBDA(LO i, LO & next, bool final) { + if (owned(i)) { + if (final) { + owned_indices(next) = i; + } + ++next; + } + }); + + out.owned_coords_2d = + Kokkos::View("owned_coords", out.num_owned, dim); + const auto owned_coords = out.owned_coords_2d; + Kokkos::View owned_gids("owned_gids", out.num_owned); + Kokkos::parallel_for( + "pcms::GatherOwnedGidsAndCoords", + Kokkos::RangePolicy(0, out.num_owned), KOKKOS_LAMBDA(LO o) { + const LO i = owned_indices(o); + owned_gids(o) = gids(i); + for (int d = 0; d < dim; ++d) { + owned_coords(o, d) = coords(i, d); + } + }); + + out.owned_indices_host = + Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out.owned_indices); + out.owned_gids_host = + Kokkos::create_mirror_view_and_copy(HostMemorySpace(), owned_gids); + return out; +} + } // namespace pcms #endif // PCMS_FIELD_LAYOUT_H diff --git a/src/pcms/field/function_space.h b/src/pcms/field/function_space.h index 78aeb8a5..67db2a08 100644 --- a/src/pcms/field/function_space.h +++ b/src/pcms/field/function_space.h @@ -151,6 +151,7 @@ inline EvaluationRequest EvaluationRequest::FromLayout( if (layout == nullptr) { throw pcms_error("EvaluationRequest::FromLayout: layout must not be null"); } + // The layout's owned DOF holders: the rows its fields' setters write. // Must evaluate GetDOFHolderCoordinates() before std::move(layout) auto coords = layout->GetDOFHolderCoordinates(); return EvaluationRequest(coords, std::move(layout), policy); diff --git a/src/pcms/field/function_space/lagrange.cpp b/src/pcms/field/function_space/lagrange.cpp index faeba79e..4fe92bdf 100644 --- a/src/pcms/field/function_space/lagrange.cpp +++ b/src/pcms/field/function_space/lagrange.cpp @@ -11,6 +11,7 @@ #include "pcms/field/evaluator/omega_h_lagrange.h" #include "pcms/field/layout/uniform_grid.h" #include "pcms/field/evaluator/uniform_grid.h" +#include "pcms/field/data/omega_h.h" #include "pcms/field/data/simple.h" #include "pcms/utility/uniform_grid.h" @@ -35,7 +36,13 @@ void ValidateLagrangeWrappedFieldData(const FieldLayout& layout, } } else #endif - { + if (dynamic_cast(&layout) != nullptr) { + if (dynamic_cast*>(&data) == nullptr) { + throw pcms_error( + "LagrangeFunctionSpace::CreateField: OmegaH layout requires " + "OmegaHFieldData"); + } + } else { if (dynamic_cast*>(&data) == nullptr) { throw pcms_error( "LagrangeFunctionSpace::CreateField: this backend requires " @@ -43,7 +50,7 @@ void ValidateLagrangeWrappedFieldData(const FieldLayout& layout, } } - if (data.GetDOFAccess()->GetOwnedValues().size() != + if (data.GetDOFAccess()->GetValues().size() != detail::ExpectedFlatFieldDataSize(layout)) { throw pcms_error( "LagrangeFunctionSpace::CreateField: field data size does not match " @@ -133,7 +140,7 @@ std::shared_ptr LagrangeFunctionSpace::FromMesh( if constexpr (std::is_same_v) { throw pcms_error("LagrangeFunctionSpace: int8_t is not supported"); } else { - return std::make_unique>(layout, basis); + return std::make_unique>(layout, basis); } }); }, @@ -170,7 +177,7 @@ std::shared_ptr LagrangeFunctionSpace::FromMesh( if constexpr (std::is_same_v) { throw pcms_error("LagrangeFunctionSpace: int8_t is not supported"); } else { - return std::make_unique>(layout, basis); + return std::make_unique>(layout, basis); } }); }, diff --git a/src/pcms/field/function_space/mfem.hpp b/src/pcms/field/function_space/mfem.hpp index a3b8820d..8a8354f9 100644 --- a/src/pcms/field/function_space/mfem.hpp +++ b/src/pcms/field/function_space/mfem.hpp @@ -74,7 +74,7 @@ class MFEMFieldFactory : public FieldFactory throw pcms_error( "MFEMFieldFactory::CreateField: requires MFEMVertexFieldData"); } - if ((*real)->GetDOFAccess()->GetOwnedValues().size() != + if ((*real)->GetDOFAccess()->GetValues().size() != detail::ExpectedFlatFieldDataSize(*layout_)) { throw pcms_error("MFEMFieldFactory::CreateField: field data size does " "not match layout"); diff --git a/src/pcms/field/function_space/polynomial_reconstruction.cpp b/src/pcms/field/function_space/polynomial_reconstruction.cpp index 2097ead7..4642c464 100644 --- a/src/pcms/field/function_space/polynomial_reconstruction.cpp +++ b/src/pcms/field/function_space/polynomial_reconstruction.cpp @@ -111,7 +111,7 @@ FieldVariant PolynomialReconstructionFunctionSpace::CreateFieldImpl( "PolynomialReconstructionFunctionSpace::CreateField: requires " "SimpleFieldData"); } - if (fd->GetDOFAccess()->GetOwnedValues().size() != + if (fd->GetDOFAccess()->GetValues().size() != detail::ExpectedFlatFieldDataSize(*layout_)) { throw pcms_error( "PolynomialReconstructionFunctionSpace::CreateField: field data size " diff --git a/src/pcms/field/function_space/spline.cpp b/src/pcms/field/function_space/spline.cpp index f169a9f4..bf98829c 100644 --- a/src/pcms/field/function_space/spline.cpp +++ b/src/pcms/field/function_space/spline.cpp @@ -61,7 +61,7 @@ FieldVariant SplineFunctionSpace::CreateFieldImpl(FieldDataVariant data) const throw pcms_error( "SplineFunctionSpace::CreateField: requires SimpleFieldData"); } - if (fd->GetDOFAccess()->GetOwnedValues().size() != + if (fd->GetDOFAccess()->GetValues().size() != detail::ExpectedFlatFieldDataSize(*layout_)) { throw pcms_error( "SplineFunctionSpace::CreateField: field data size does not match " diff --git a/src/pcms/field/function_space/xgc.h b/src/pcms/field/function_space/xgc.h index ee0b63b2..94198e64 100644 --- a/src/pcms/field/function_space/xgc.h +++ b/src/pcms/field/function_space/xgc.h @@ -63,7 +63,7 @@ class XGCFieldFactory : public FieldFactory throw pcms_error( "XGCFieldFactory::CreateField: requires XGCFieldData"); } - if (fd->GetDOFAccess()->GetOwnedValues().size() != + if (fd->GetDOFAccess()->GetValues().size() != static_cast(layout_->GetFullDataSize())) { throw pcms_error( "XGCFieldFactory::CreateField: field data size does not match " diff --git a/src/pcms/field/layout/empty.cpp b/src/pcms/field/layout/empty.cpp index cdd6dc0b..fe5ad996 100644 --- a/src/pcms/field/layout/empty.cpp +++ b/src/pcms/field/layout/empty.cpp @@ -56,7 +56,7 @@ int EmptyFieldLayout::GetNumComponents() const return num_components_; } -LO EmptyFieldLayout::GetNumOwnedDofHolder() const +LO EmptyFieldLayout::GetNumDofHolder() const { return 0; } diff --git a/src/pcms/field/layout/empty.h b/src/pcms/field/layout/empty.h index c3e12970..85031795 100644 --- a/src/pcms/field/layout/empty.h +++ b/src/pcms/field/layout/empty.h @@ -21,7 +21,7 @@ class EmptyFieldLayout : public FieldLayout const noexcept override; int GetNumComponents() const override; - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; GlobalIDView GetGidsHost() const override; diff --git a/src/pcms/field/layout/mesh_fields.cpp b/src/pcms/field/layout/mesh_fields.cpp index e93e02bc..caeaef56 100644 --- a/src/pcms/field/layout/mesh_fields.cpp +++ b/src/pcms/field/layout/mesh_fields.cpp @@ -208,7 +208,7 @@ MeshFieldsAdapterLayout::MeshFieldsAdapterLayout( global_id_name_(global_id_name), num_components_(num_components), nodes_per_dim_(nodes_per_dim), - dof_holder_coords_("", GetNumOwnedDofHolder(), mesh_.dim()), + dof_holder_coords_("", GetNumDofHolderWithGhosts(), mesh_.dim()), class_ids_(GetNumEnts()), class_dims_(class_ids_.size()), owned_("", class_dims_.size()), @@ -246,21 +246,10 @@ MeshFieldsAdapterLayout::MeshFieldsAdapterLayout( offset += mesh.nents(i); } } - gids_host_ = Omega_h::HostWrite(gids_); - - int n = class_ids_.size(); - classification_dims_host_ = - Kokkos::View("classification_dims", n); - classification_ids_host_ = - Kokkos::View("classification_ids", n); - auto class_dims_h = Omega_h::HostRead(class_dims_); - auto class_ids_h = Omega_h::HostRead(class_ids_); - for (int i = 0; i < n; ++i) { - classification_dims_host_(i) = - static_cast(static_cast(class_dims_h[i])); - classification_ids_host_(i) = static_cast(class_ids_h[i]); - } + BuildOwnedViews(); + discretization_ = std::make_shared(mesh_); + BuildOwnedClassification(); } std::shared_ptr @@ -274,7 +263,7 @@ int MeshFieldsAdapterLayout::GetNumComponents() const return num_components_; } -LO MeshFieldsAdapterLayout::GetNumOwnedDofHolder() const +LO MeshFieldsAdapterLayout::GetNumDofHolderWithGhosts() const { LO count = 0; for (int i = 0; i <= mesh_.dim(); ++i) { @@ -283,6 +272,16 @@ LO MeshFieldsAdapterLayout::GetNumOwnedDofHolder() const return count; } +LO MeshFieldsAdapterLayout::GetNumDofHolder() const +{ + return num_owned_; +} + +bool MeshFieldsAdapterLayout::HasGhosts() const +{ + return num_owned_ != GetNumDofHolderWithGhosts(); +} + GO MeshFieldsAdapterLayout::GetNumGlobalDofHolder() const { LO count = 0; @@ -306,16 +305,44 @@ Rank1View MeshFieldsAdapterLayout::GetOwnedHost() GlobalIDView MeshFieldsAdapterLayout::GetGidsHost() const { - return GlobalIDView(gids_host_.data(), gids_host_.size()); + return GlobalIDView(owned_gids_host_.data(), + owned_gids_host_.size()); } CoordinateView MeshFieldsAdapterLayout::GetDOFHolderCoordinates() const { - auto coords_view = MakeConstRank2View(dof_holder_coords_); + auto coords_view = MakeConstRank2View(owned_coords_2d_); return CoordinateView{GetCoordinateSystem(), coords_view}; } +Kokkos::View +MeshFieldsAdapterLayout::GetOwnedIndicesHost() const +{ + return owned_indices_host_; +} + +Kokkos::View +MeshFieldsAdapterLayout::GetOwnedIndices() const +{ + return owned_indices_; +} + +void MeshFieldsAdapterLayout::BuildOwnedViews() +{ + Kokkos::deep_copy(owned_host_, owned_); + + auto owned = BuildOwnedLayoutData( + owned_, + Kokkos::View(gids_.data(), gids_.size()), + dof_holder_coords_, mesh_.dim()); + num_owned_ = owned.num_owned; + owned_indices_host_ = owned.owned_indices_host; + owned_gids_host_ = owned.owned_gids_host; + owned_coords_2d_ = owned.owned_coords_2d; + owned_indices_ = owned.owned_indices; +} + bool MeshFieldsAdapterLayout::IsDistributed() const { return true; @@ -349,6 +376,13 @@ Omega_h::Mesh& MeshFieldsAdapterLayout::GetMesh() const } EntOffsetsArray MeshFieldsAdapterLayout::GetEntOffsets() const +{ + return !HasGhosts() + ? EntOffsetsWithGhosts() + : OwnedEntOffsets(EntOffsetsWithGhosts(), owned_indices_host_); +} + +EntOffsetsArray MeshFieldsAdapterLayout::EntOffsetsWithGhosts() const { EntOffsetsArray offsets{}; size_t offset = 0; @@ -366,16 +400,39 @@ int MeshFieldsAdapterLayout::GetDimension() const return mesh_.dim(); } +int MeshFieldsAdapterLayout::GetDOFHolderEntityDim() const +{ + // One DOF holder per entity of a single dimension, or none. + int entity_dim = -1; + for (int i = 0; i <= mesh_.dim(); ++i) { + if (nodes_per_dim_[i] == 0) { + continue; + } + if (entity_dim >= 0 || nodes_per_dim_[i] != 1) { + return -1; + } + entity_dim = i; + } + return entity_dim; +} + Rank1View MeshFieldsAdapterLayout::GetDOFHolderClassificationDimensionsHost() const { - return make_const_array_view(classification_dims_host_); + return make_const_array_view(owned_classification_dims_host_); } Rank1View MeshFieldsAdapterLayout::GetDOFHolderClassificationIdsHost() const { - return make_const_array_view(classification_ids_host_); + return make_const_array_view(owned_classification_ids_host_); +} + +void MeshFieldsAdapterLayout::BuildOwnedClassification() +{ + BuildOwnedClassificationHost(class_dims_, class_ids_, owned_indices_, + owned_classification_dims_host_, + owned_classification_ids_host_); } } // namespace pcms diff --git a/src/pcms/field/layout/mesh_fields.h b/src/pcms/field/layout/mesh_fields.h index 004fd202..1c7f2586 100644 --- a/src/pcms/field/layout/mesh_fields.h +++ b/src/pcms/field/layout/mesh_fields.h @@ -25,24 +25,29 @@ class MeshFieldsAdapterLayout : public FieldLayout int GetNumComponents() const override; // nodes for standard lagrange FEM - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; + [[nodiscard]] bool HasGhosts() const override; + // Number of DOF holders, owned and ghost: the size of the field data + // this layout's backends store. + LO GetNumDofHolderWithGhosts() const; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; GlobalIDView GetGidsHost() const override; CoordinateView GetDOFHolderCoordinates() const override; + EntOffsetsArray GetEntOffsets() const override; + Kokkos::View GetOwnedIndicesHost() const override; + Kokkos::View GetOwnedIndices() const override; - // returns true if the field layout is distributed - // if the field layout is distributed, the owned and global dofs are the same + // returns true if the field layout is distributed: on several ranks it may + // hold ghost DOF holders in addition to the owned ones (see HasGhosts) [[nodiscard]] bool IsDistributed() const override; - EntOffsetsArray GetEntOffsets() const override; - int GetDimension() const override; + int GetDOFHolderEntityDim() const override; Rank1View GetDOFHolderClassificationDimensionsHost() const override; - Rank1View GetDOFHolderClassificationIdsHost() const override; @@ -56,7 +61,6 @@ class MeshFieldsAdapterLayout : public FieldLayout Omega_h::Mesh& mesh_; Omega_h::Write gids_; - Omega_h::HostWrite gids_host_; std::string global_id_name_; int num_components_; std::array nodes_per_dim_; @@ -67,8 +71,21 @@ class MeshFieldsAdapterLayout : public FieldLayout Omega_h::HostWrite class_dims_host_; Kokkos::View owned_; Kokkos::View owned_host_; - Kokkos::View classification_dims_host_; - Kokkos::View classification_ids_host_; + LO num_owned_ = 0; + Kokkos::View owned_gids_host_; + Kokkos::View owned_coords_2d_; + Kokkos::View owned_indices_host_; + Kokkos::View owned_indices_; + + Kokkos::View owned_classification_dims_host_; + Kokkos::View owned_classification_ids_host_; + + // Entity offsets over the DOF holders with ghosts. + EntOffsetsArray EntOffsetsWithGhosts() const; + void BuildOwnedViews(); + // Requires the owned views and the classification arrays with ghosts. + void BuildOwnedClassification(); + std::shared_ptr discretization_; }; diff --git a/src/pcms/field/layout/mfem.cpp b/src/pcms/field/layout/mfem.cpp index 3e70e584..d98e3574 100644 --- a/src/pcms/field/layout/mfem.cpp +++ b/src/pcms/field/layout/mfem.cpp @@ -86,7 +86,7 @@ int MFEMLayout::GetNumComponents() const return 1; } -LO MFEMLayout::GetNumOwnedDofHolder() const +LO MFEMLayout::GetNumDofHolder() const { return static_cast(pmesh_.GetNV()); } @@ -121,7 +121,7 @@ EntOffsetsArray MFEMLayout::GetEntOffsets() const { EntOffsetsArray offsets{}; offsets.fill(0); - const auto n = static_cast(GetNumOwnedDofHolder()); + const auto n = static_cast(GetNumDofHolder()); for (int i = 1; i < ent_offsets_len; ++i) { offsets[i] = n; } diff --git a/src/pcms/field/layout/mfem.hpp b/src/pcms/field/layout/mfem.hpp index edad78ea..de0e3585 100644 --- a/src/pcms/field/layout/mfem.hpp +++ b/src/pcms/field/layout/mfem.hpp @@ -22,7 +22,7 @@ class MFEMLayout : public FieldLayout const noexcept override; int GetNumComponents() const override; - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; diff --git a/src/pcms/field/layout/omega_h_entity.cpp b/src/pcms/field/layout/omega_h_entity.cpp index bbf5e4ae..7aee5924 100644 --- a/src/pcms/field/layout/omega_h_entity.cpp +++ b/src/pcms/field/layout/omega_h_entity.cpp @@ -67,26 +67,18 @@ OmegaHEntityLayout::OmegaHEntityLayout( class_dims_(mesh.get_array(entity_dim, "class_dim")), owned_(BuildOwned(mesh, entity_dim)), owned_host_("owned_host", owned_.size()), - classification_dims_host_("classification_dims", mesh.nents(entity_dim)), - classification_ids_host_("classification_ids", mesh.nents(entity_dim)), discretization_(std::make_shared(mesh)) { SetCoordinateSystem( ResolveCoordinateSystem(std::move(coordinate_system), dimension_)); PCMS_ALWAYS_ASSERT(entity_dim_ >= 0 && entity_dim_ <= dimension_); - gids_host_ = Omega_h::HostWrite(gids_); Kokkos::deep_copy(owned_host_, owned_); + BuildOwnedViews(); class_dims_ = Omega_h::Read(class_dims_); class_ids_ = Omega_h::Read(class_ids_); - auto class_dims_host = Omega_h::HostRead(class_dims_); - auto class_ids_host = Omega_h::HostRead(class_ids_); - for (int i = 0; i < mesh.nents(entity_dim_); ++i) { - classification_dims_host_(i) = - static_cast(static_cast(class_dims_host[i])); - classification_ids_host_(i) = static_cast(class_ids_host[i]); - } + BuildOwnedClassification(); } std::shared_ptr OmegaHEntityLayout::GetDiscretization() @@ -100,11 +92,21 @@ int OmegaHEntityLayout::GetNumComponents() const return num_components_; } -LO OmegaHEntityLayout::GetNumOwnedDofHolder() const +LO OmegaHEntityLayout::GetNumDofHolderWithGhosts() const { return static_cast(coords_.size() / dimension_); } +LO OmegaHEntityLayout::GetNumDofHolder() const +{ + return num_owned_; +} + +bool OmegaHEntityLayout::HasGhosts() const +{ + return num_owned_ != GetNumDofHolderWithGhosts(); +} + GO OmegaHEntityLayout::GetNumGlobalDofHolder() const { return num_global_dof_holder_; @@ -117,14 +119,44 @@ Rank1View OmegaHEntityLayout::GetOwnedHost() const GlobalIDView OmegaHEntityLayout::GetGidsHost() const { - return GlobalIDView(gids_host_.data(), gids_host_.size()); + return GlobalIDView(owned_gids_host_.data(), + owned_gids_host_.size()); } CoordinateView OmegaHEntityLayout::GetDOFHolderCoordinates() const { - return CoordinateView{GetCoordinateSystem(), - MakeConstRank2View(coords_2d_)}; + using LayoutPolicy = + detail::default_layout_for_memory_space_t; + Rank2View coords_view( + owned_coords_2d_.data(), num_owned_, dimension_); + return CoordinateView{GetCoordinateSystem(), + coords_view}; +} + +Kokkos::View +OmegaHEntityLayout::GetOwnedIndicesHost() const +{ + return owned_indices_host_; +} + +Kokkos::View OmegaHEntityLayout::GetOwnedIndices() + const +{ + return owned_indices_; +} + +void OmegaHEntityLayout::BuildOwnedViews() +{ + auto owned = BuildOwnedLayoutData( + owned_, + Kokkos::View(gids_.data(), gids_.size()), + coords_2d_, dimension_); + num_owned_ = owned.num_owned; + owned_indices_host_ = owned.owned_indices_host; + owned_gids_host_ = owned.owned_gids_host; + owned_coords_2d_ = owned.owned_coords_2d; + owned_indices_ = owned.owned_indices; } bool OmegaHEntityLayout::IsDistributed() const @@ -133,10 +165,17 @@ bool OmegaHEntityLayout::IsDistributed() const } EntOffsetsArray OmegaHEntityLayout::GetEntOffsets() const +{ + return !HasGhosts() + ? EntOffsetsWithGhosts() + : OwnedEntOffsets(EntOffsetsWithGhosts(), owned_indices_host_); +} + +EntOffsetsArray OmegaHEntityLayout::EntOffsetsWithGhosts() const { EntOffsetsArray offsets{}; offsets.fill(0); - const auto n = static_cast(GetNumOwnedDofHolder()); + const auto n = static_cast(GetNumDofHolderWithGhosts()); for (int i = entity_dim_ + 1; i < ent_offsets_len; ++i) offsets[i] = n; return offsets; @@ -147,16 +186,28 @@ int OmegaHEntityLayout::GetDimension() const return dimension_; } +int OmegaHEntityLayout::GetDOFHolderEntityDim() const +{ + return entity_dim_; +} + Rank1View OmegaHEntityLayout::GetDOFHolderClassificationDimensionsHost() const { - return make_const_array_view(classification_dims_host_); + return make_const_array_view(owned_classification_dims_host_); } Rank1View OmegaHEntityLayout::GetDOFHolderClassificationIdsHost() const { - return make_const_array_view(classification_ids_host_); + return make_const_array_view(owned_classification_ids_host_); +} + +void OmegaHEntityLayout::BuildOwnedClassification() +{ + BuildOwnedClassificationHost(class_dims_, class_ids_, owned_indices_, + owned_classification_dims_host_, + owned_classification_ids_host_); } } // namespace pcms diff --git a/src/pcms/field/layout/omega_h_entity.h b/src/pcms/field/layout/omega_h_entity.h index 65f65147..dba0f84a 100644 --- a/src/pcms/field/layout/omega_h_entity.h +++ b/src/pcms/field/layout/omega_h_entity.h @@ -20,20 +20,26 @@ class OmegaHEntityLayout : public FieldLayout const noexcept override; int GetNumComponents() const override; - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; + [[nodiscard]] bool HasGhosts() const override; + // Number of DOF holders, owned and ghost: the size of the field data + // this layout's backends store. + LO GetNumDofHolderWithGhosts() const; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; GlobalIDView GetGidsHost() const override; CoordinateView GetDOFHolderCoordinates() const override; + EntOffsetsArray GetEntOffsets() const override; + Kokkos::View GetOwnedIndicesHost() const override; + Kokkos::View GetOwnedIndices() const override; [[nodiscard]] bool IsDistributed() const override; - EntOffsetsArray GetEntOffsets() const override; int GetDimension() const override; + int GetDOFHolderEntityDim() const override; Rank1View GetDOFHolderClassificationDimensionsHost() const override; - Rank1View GetDOFHolderClassificationIdsHost() const override; @@ -44,15 +50,27 @@ class OmegaHEntityLayout : public FieldLayout GO num_global_dof_holder_; Omega_h::Write gids_; - Omega_h::HostWrite gids_host_; Omega_h::Read coords_; // device coordinates (1D flattened) Kokkos::View coords_2d_; // device coordinates (2D) Omega_h::Read class_ids_; Omega_h::Read class_dims_; Kokkos::View owned_; Kokkos::View owned_host_; - Kokkos::View classification_dims_host_; - Kokkos::View classification_ids_host_; + LO num_owned_ = 0; + Kokkos::View owned_gids_host_; + Kokkos::View owned_coords_2d_; + Kokkos::View owned_indices_host_; + Kokkos::View owned_indices_; + + Kokkos::View owned_classification_dims_host_; + Kokkos::View owned_classification_ids_host_; + + // Entity offsets over the DOF holders with ghosts. + EntOffsetsArray EntOffsetsWithGhosts() const; + void BuildOwnedViews(); + // Requires the owned views and the classification arrays with ghosts. + void BuildOwnedClassification(); + std::shared_ptr discretization_; }; diff --git a/src/pcms/field/layout/omega_h_lagrange.cpp b/src/pcms/field/layout/omega_h_lagrange.cpp index df1b6f25..f0902bbe 100644 --- a/src/pcms/field/layout/omega_h_lagrange.cpp +++ b/src/pcms/field/layout/omega_h_lagrange.cpp @@ -90,35 +90,25 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout( EnsureClassification(mesh_); gids_ = BuildGids(mesh_, entity_dim, global_id_name_); - gids_host_ = Omega_h::HostWrite(gids_); coords_ = get_entity_centroids(mesh_, entity_dim); coords_2d_ = ConvertCoordsTo2D(coords_, mesh_.nents(entity_dim), mesh_.dim()); owned_ = BuildOwned(mesh_, entity_dim); owned_host_ = Kokkos::View("owned_host", owned_.size()); Kokkos::deep_copy(owned_host_, owned_); + BuildOwnedViews(); class_ids_ = Omega_h::Read( mesh_.get_array(entity_dim, "class_id")); class_dims_ = Omega_h::Read( mesh_.get_array(entity_dim, "class_dim")); - auto class_ids_host = Omega_h::HostRead(class_ids_); - auto class_dims_host = Omega_h::HostRead(class_dims_); - - int n = mesh_.nents(entity_dim); - classification_dims_host_ = - Kokkos::View("classification_dims", n); - classification_ids_host_ = - Kokkos::View("classification_ids", n); - for (int i = 0; i < n; ++i) { - classification_dims_host_(i) = - static_cast(static_cast(class_dims_host[i])); - classification_ids_host_(i) = static_cast(class_ids_host[i]); - } discretization_ = std::make_shared(mesh_); + BuildOwnedClassification(); - BuildGlobalToLocalPermutation(); + BuildGlobalToLocalPermutation( + Kokkos::View(gids_.data(), gids_.size()), + EntOffsetsWithGhosts()); } OmegaHLagrangeLayout::OmegaHLagrangeLayout( @@ -137,35 +127,25 @@ OmegaHLagrangeLayout::OmegaHLagrangeLayout( EnsureClassification(mesh_); gids_ = BuildGids(mesh_, entity_dim, global_id_name_); - gids_host_ = Omega_h::HostWrite(gids_); coords_ = get_entity_centroids(mesh_, entity_dim); coords_2d_ = ConvertCoordsTo2D(coords_, mesh_.nents(entity_dim), mesh_.dim()); owned_ = BuildOwned(mesh_, entity_dim, owned_mask); owned_host_ = Kokkos::View("owned_host", owned_.size()); Kokkos::deep_copy(owned_host_, owned_); + BuildOwnedViews(); class_ids_ = Omega_h::Read( mesh_.get_array(entity_dim, "class_id")); class_dims_ = Omega_h::Read( mesh_.get_array(entity_dim, "class_dim")); - auto class_ids_host = Omega_h::HostRead(class_ids_); - auto class_dims_host = Omega_h::HostRead(class_dims_); - - int n = mesh_.nents(entity_dim); - classification_dims_host_ = - Kokkos::View("classification_dims", n); - classification_ids_host_ = - Kokkos::View("classification_ids", n); - for (int i = 0; i < n; ++i) { - classification_dims_host_(i) = - static_cast(static_cast(class_dims_host[i])); - classification_ids_host_(i) = static_cast(class_ids_host[i]); - } discretization_ = std::make_shared(mesh_); + BuildOwnedClassification(); - BuildGlobalToLocalPermutation(); + BuildGlobalToLocalPermutation( + Kokkos::View(gids_.data(), gids_.size()), + EntOffsetsWithGhosts()); } std::shared_ptr OmegaHLagrangeLayout::GetDiscretization() @@ -179,11 +159,21 @@ int OmegaHLagrangeLayout::GetNumComponents() const return num_components_; } -LO OmegaHLagrangeLayout::GetNumOwnedDofHolder() const +LO OmegaHLagrangeLayout::GetNumDofHolderWithGhosts() const { return mesh_.nents(EntityDimForOrder(order_, mesh_.dim())); } +LO OmegaHLagrangeLayout::GetNumDofHolder() const +{ + return num_owned_; +} + +bool OmegaHLagrangeLayout::HasGhosts() const +{ + return num_owned_ != GetNumDofHolderWithGhosts(); +} + GO OmegaHLagrangeLayout::GetNumGlobalDofHolder() const { return mesh_.nglobal_ents(EntityDimForOrder(order_, mesh_.dim())); @@ -197,7 +187,8 @@ Rank1View OmegaHLagrangeLayout::GetOwnedHost() GlobalIDView OmegaHLagrangeLayout::GetGidsHost() const { - return GlobalIDView(gids_host_.data(), gids_host_.size()); + return GlobalIDView(owned_gids_host_.data(), + owned_gids_host_.size()); } GlobalIDView OmegaHLagrangeLayout::GetGids() const @@ -208,8 +199,37 @@ GlobalIDView OmegaHLagrangeLayout::GetGids() const CoordinateView OmegaHLagrangeLayout::GetDOFHolderCoordinates() const { - return CoordinateView{GetCoordinateSystem(), - MakeConstRank2View(coords_2d_)}; + using LayoutPolicy = + detail::default_layout_for_memory_space_t; + Rank2View coords_view( + owned_coords_2d_.data(), num_owned_, mesh_.dim()); + return CoordinateView{GetCoordinateSystem(), + coords_view}; +} + +Kokkos::View +OmegaHLagrangeLayout::GetOwnedIndicesHost() const +{ + return owned_indices_host_; +} + +Kokkos::View +OmegaHLagrangeLayout::GetOwnedIndices() const +{ + return owned_indices_; +} + +void OmegaHLagrangeLayout::BuildOwnedViews() +{ + auto owned = BuildOwnedLayoutData( + owned_, + Kokkos::View(gids_.data(), gids_.size()), + coords_2d_, mesh_.dim()); + num_owned_ = owned.num_owned; + owned_indices_host_ = owned.owned_indices_host; + owned_gids_host_ = owned.owned_gids_host; + owned_coords_2d_ = owned.owned_coords_2d; + owned_indices_ = owned.owned_indices; } bool OmegaHLagrangeLayout::IsDistributed() const @@ -218,6 +238,13 @@ bool OmegaHLagrangeLayout::IsDistributed() const } EntOffsetsArray OmegaHLagrangeLayout::GetEntOffsets() const +{ + return !HasGhosts() + ? EntOffsetsWithGhosts() + : OwnedEntOffsets(EntOffsetsWithGhosts(), owned_indices_host_); +} + +EntOffsetsArray OmegaHLagrangeLayout::EntOffsetsWithGhosts() const { // Slot i holds the starting DOF index for entity dimension i. // Slot 4 holds the total DOF count. @@ -237,16 +264,28 @@ int OmegaHLagrangeLayout::GetDimension() const return mesh_.dim(); } +int OmegaHLagrangeLayout::GetDOFHolderEntityDim() const +{ + return EntityDimForOrder(order_, mesh_.dim()); +} + Rank1View OmegaHLagrangeLayout::GetDOFHolderClassificationDimensionsHost() const { - return make_const_array_view(classification_dims_host_); + return make_const_array_view(owned_classification_dims_host_); } Rank1View OmegaHLagrangeLayout::GetDOFHolderClassificationIdsHost() const { - return make_const_array_view(classification_ids_host_); + return make_const_array_view(owned_classification_ids_host_); +} + +void OmegaHLagrangeLayout::BuildOwnedClassification() +{ + BuildOwnedClassificationHost(class_dims_, class_ids_, owned_indices_, + owned_classification_dims_host_, + owned_classification_ids_host_); } int OmegaHLagrangeLayout::GetOrder() const diff --git a/src/pcms/field/layout/omega_h_lagrange.h b/src/pcms/field/layout/omega_h_lagrange.h index 7951bb0e..98e354dc 100644 --- a/src/pcms/field/layout/omega_h_lagrange.h +++ b/src/pcms/field/layout/omega_h_lagrange.h @@ -31,23 +31,28 @@ class OmegaHLagrangeLayout : public FieldLayout const noexcept override; int GetNumComponents() const override; - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; + [[nodiscard]] bool HasGhosts() const override; + // Number of DOF holders, owned and ghost: the size of the field data + // this layout's backends store. + LO GetNumDofHolderWithGhosts() const; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; GlobalIDView GetGidsHost() const override; GlobalIDView GetGids() const; CoordinateView GetDOFHolderCoordinates() const override; + EntOffsetsArray GetEntOffsets() const override; + Kokkos::View GetOwnedIndicesHost() const override; + Kokkos::View GetOwnedIndices() const override; [[nodiscard]] bool IsDistributed() const override; - EntOffsetsArray GetEntOffsets() const override; - int GetDimension() const override; + int GetDOFHolderEntityDim() const override; Rank1View GetDOFHolderClassificationDimensionsHost() const override; - Rank1View GetDOFHolderClassificationIdsHost() const override; @@ -61,15 +66,27 @@ class OmegaHLagrangeLayout : public FieldLayout std::string global_id_name_; Omega_h::Write gids_; - Omega_h::HostWrite gids_host_; Kokkos::View coords_2d_; Omega_h::Read coords_; // device coordinates Kokkos::View owned_; Kokkos::View owned_host_; Omega_h::Read class_ids_; Omega_h::Read class_dims_; - Kokkos::View classification_dims_host_; - Kokkos::View classification_ids_host_; + LO num_owned_ = 0; + Kokkos::View owned_gids_host_; + Kokkos::View owned_coords_2d_; + Kokkos::View owned_indices_host_; + Kokkos::View owned_indices_; + + Kokkos::View owned_classification_dims_host_; + Kokkos::View owned_classification_ids_host_; + + // Entity offsets over the DOF holders with ghosts. + EntOffsetsArray EntOffsetsWithGhosts() const; + void BuildOwnedViews(); + // Requires the owned views and the classification arrays. + void BuildOwnedClassification(); + std::shared_ptr discretization_; }; diff --git a/src/pcms/field/layout/point_cloud.cpp b/src/pcms/field/layout/point_cloud.cpp index 1c020e27..51a7aba9 100644 --- a/src/pcms/field/layout/point_cloud.cpp +++ b/src/pcms/field/layout/point_cloud.cpp @@ -95,7 +95,7 @@ int PointCloudLayout::GetNumComponents() const return components_; } -LO PointCloudLayout::GetNumOwnedDofHolder() const +LO PointCloudLayout::GetNumDofHolder() const { return coords_.extent(0); } diff --git a/src/pcms/field/layout/point_cloud.h b/src/pcms/field/layout/point_cloud.h index 1754c0ad..29d2d0bd 100644 --- a/src/pcms/field/layout/point_cloud.h +++ b/src/pcms/field/layout/point_cloud.h @@ -22,7 +22,7 @@ class PointCloudLayout : public FieldLayout int GetNumComponents() const override; // nodes for standard lagrange FEM - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; diff --git a/src/pcms/field/layout/uniform_grid.cpp b/src/pcms/field/layout/uniform_grid.cpp index 5887d096..926af7df 100644 --- a/src/pcms/field/layout/uniform_grid.cpp +++ b/src/pcms/field/layout/uniform_grid.cpp @@ -125,18 +125,18 @@ UniformGridFieldLayout::UniformGridFieldLayout( : grid_(std::move(grid)), num_components_(num_components), order_(order), - gids_("gids", GetNumDofHolders()), - gids_host_("gids_host", GetNumDofHolders()), - dof_holder_coords_("dof_holder_coords", GetNumDofHolders(), Dim), - owned_("owned", GetNumDofHolders()), - owned_host_("owned_host", GetNumDofHolders()) + gids_("gids", CountDofHolders()), + gids_host_("gids_host", CountDofHolders()), + dof_holder_coords_("dof_holder_coords", CountDofHolders(), Dim), + owned_("owned", CountDofHolders()), + owned_host_("owned_host", CountDofHolders()) { PCMS_FUNCTION_TIMER; SetCoordinateSystem(ResolveCoordinateSystem(std::move(coordinate_system), static_cast(Dim))); PCMS_ALWAYS_ASSERT(order_ == 0 || order_ == 1); - LO num_dofs = GetNumDofHolders(); + LO num_dofs = CountDofHolders(); Kokkos::parallel_for( "InitUniformGridGidsAndOwned", @@ -175,7 +175,7 @@ UniformGridFieldLayout::UniformGridFieldLayout( } int entity_dim = (order_ == 0) ? static_cast(Dim) : 0; - LO n = GetNumDofHolders(); + LO n = CountDofHolders(); classification_dims_ = Kokkos::View("classification_dims", n); classification_ids_ = @@ -206,15 +206,15 @@ int UniformGridFieldLayout::GetNumComponents() const } template -LO UniformGridFieldLayout::GetNumOwnedDofHolder() const +LO UniformGridFieldLayout::GetNumDofHolder() const { - return GetNumDofHolders(); + return CountDofHolders(); } template GO UniformGridFieldLayout::GetNumGlobalDofHolder() const { - return GetNumDofHolders(); + return CountDofHolders(); } template @@ -267,7 +267,7 @@ LO UniformGridFieldLayout::GetNumVertices() const } template -LO UniformGridFieldLayout::GetNumDofHolders() const +LO UniformGridFieldLayout::CountDofHolders() const { return order_ == 0 ? GetNumCells() : GetNumVertices(); } @@ -297,7 +297,7 @@ EntOffsetsArray UniformGridFieldLayout::GetEntOffsets() const { EntOffsetsArray offsets{}; offsets.fill(0); - LO n = GetNumDofHolders(); + LO n = CountDofHolders(); int entity_dim = (order_ == 0) ? static_cast(Dim) : 0; for (int i = entity_dim + 1; i < ent_offsets_len; ++i) { offsets[i] = n; diff --git a/src/pcms/field/layout/uniform_grid.h b/src/pcms/field/layout/uniform_grid.h index 21828edd..ff0cf561 100644 --- a/src/pcms/field/layout/uniform_grid.h +++ b/src/pcms/field/layout/uniform_grid.h @@ -23,7 +23,7 @@ class UniformGridFieldLayout : public FieldLayout const noexcept override; int GetNumComponents() const override; - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; @@ -48,7 +48,7 @@ class UniformGridFieldLayout : public FieldLayout int GetOrder() const; private: - LO GetNumDofHolders() const; + LO CountDofHolders() const; UniformGrid grid_; int num_components_; diff --git a/src/pcms/field/layout/xgc.cpp b/src/pcms/field/layout/xgc.cpp index 54d9eecf..8f78a6f5 100644 --- a/src/pcms/field/layout/xgc.cpp +++ b/src/pcms/field/layout/xgc.cpp @@ -89,7 +89,7 @@ int XGCFieldLayout::GetNumComponents() const return 1; } -LO XGCFieldLayout::GetNumOwnedDofHolder() const +LO XGCFieldLayout::GetNumDofHolder() const { return num_plane_nodes_; } diff --git a/src/pcms/field/layout/xgc.h b/src/pcms/field/layout/xgc.h index 8d3e3ff4..d51e072c 100644 --- a/src/pcms/field/layout/xgc.h +++ b/src/pcms/field/layout/xgc.h @@ -18,7 +18,7 @@ class XGCFieldLayout : public FieldLayout const noexcept override; int GetNumComponents() const override; - LO GetNumOwnedDofHolder() const override; + LO GetNumDofHolder() const override; GO GetNumGlobalDofHolder() const override; Rank1View GetOwnedHost() const override; GlobalIDView GetGidsHost() const override; diff --git a/src/pcms/field/uniform_grid_binary_field.h b/src/pcms/field/uniform_grid_binary_field.h index 1643815e..c02bdea6 100644 --- a/src/pcms/field/uniform_grid_binary_field.h +++ b/src/pcms/field/uniform_grid_binary_field.h @@ -46,7 +46,7 @@ CreateUniformGridBinaryField(Omega_h::Mesh& mesh, const UniformGrid& grid) auto coord_view = layout->GetDOFHolderCoordinates(); auto coords = coord_view.GetValues(); - LO n = layout->GetNumOwnedDofHolder(); + LO n = layout->GetNumDofHolder(); Kokkos::View coords_d("coords_d", n); Kokkos::parallel_for( diff --git a/src/pcms/localization/localization_path_selection.h b/src/pcms/localization/localization_path_selection.h index e08069cb..b3c0ce20 100644 --- a/src/pcms/localization/localization_path_selection.h +++ b/src/pcms/localization/localization_path_selection.h @@ -17,9 +17,10 @@ enum class LocalizationPath }; // Returns the unique entity dimension represented by all DOF holders in the -// layout when its offsets contain exactly one populated entity block and the -// full DOF-holder coordinate set matches that discretization entity count. -// Otherwise, returns std::nullopt. +// layout: the layout's GetDOFHolderEntityDim() when it has one DOF holder per +// local mesh entity of a dimension, or else, for a layout without ghosts, the +// dimension of its only populated entity block when its DOF-holder count +// matches that discretization entity count. Otherwise, returns std::nullopt. inline std::optional GetUniformDofHolderEntityDim( const FieldLayout& layout) { @@ -28,6 +29,12 @@ inline std::optional GetUniformDofHolderEntityDim( return std::nullopt; } + // A layout with one DOF holder per local mesh entity of a dimension. + if (const int mesh_entity_dim = layout.GetDOFHolderEntityDim(); + mesh_entity_dim >= 0) { + return mesh_entity_dim; + } + const auto offsets = layout.GetEntOffsets(); std::optional entity_dim; for (int dim = 0; dim < ent_offsets_len - 1; ++dim) { @@ -44,8 +51,9 @@ inline std::optional GetUniformDofHolderEntityDim( return std::nullopt; } - if (static_cast(layout.GetDOFHolderCoordinates().GetValues().extent(0)) != - disc->GetNumEntities(*entity_dim)) { + // Only a layout without ghosts covers exactly its discretization's entities. + if (layout.HasGhosts() || + layout.GetNumDofHolder() != disc->GetNumEntities(*entity_dim)) { return std::nullopt; } @@ -71,9 +79,12 @@ inline LocalizationPath SelectLocalizationPath( auto target_disc = target_layout->GetDiscretization(); auto target_entity_dim = GetUniformDofHolderEntityDim(*target_layout); - if (source_entity_dim.has_value() && *source_entity_dim == Face && - target_entity_dim.has_value() && *target_entity_dim == Vertex && - target_disc != nullptr && source_disc->GetDimension() == 2 && + // The same-mesh shortcut answers for every local target vertex, so it + // applies only when those are exactly the owned query points. + if (!target_layout->HasGhosts() && source_entity_dim.has_value() && + *source_entity_dim == Face && target_entity_dim.has_value() && + *target_entity_dim == Vertex && target_disc != nullptr && + source_disc->GetDimension() == 2 && source_disc->SameEntities(*target_disc)) { return LocalizationPath::CentroidToVertexAdjacencySearch; } diff --git a/src/pcms/pythonapi/bind_field_base.cpp b/src/pcms/pythonapi/bind_field_base.cpp index b38c8fbf..df3906e8 100644 --- a/src/pcms/pythonapi/bind_field_base.cpp +++ b/src/pcms/pythonapi/bind_field_base.cpp @@ -59,7 +59,8 @@ void bind_create_field_module(py::module& m) py::class_>( m, "FieldLayout", "How a field's DOF holders are laid out") .def("get_num_components", &FieldLayout::GetNumComponents) - .def("get_num_owned_dof_holder", &FieldLayout::GetNumOwnedDofHolder) + .def("get_num_dof_holder", &FieldLayout::GetNumDofHolder, + "Number of DOF holders this rank owns") .def("get_num_global_dof_holder", &FieldLayout::GetNumGlobalDofHolder) .def("get_dimension", &FieldLayout::GetDimension) .def("is_distributed", &FieldLayout::IsDistributed) @@ -75,7 +76,7 @@ void bind_create_field_module(py::module& m) return ToPyCoordinateView( self.cast().GetDOFHolderCoordinates(), self); }, - "The DOF-holder coordinates (num_dof_holders × dim), read-only"); + "The owned DOF-holder coordinates (num_dof_holders × dim), read-only"); py::class_>(m, "FieldData", "A field's stored values") .def("get_value_type", &FieldData::GetValueType) @@ -121,7 +122,7 @@ void bind_create_field_module(py::module& m) "set_dof_holder_data_unchecked", [](Field& self, const NumpyInput& data) { const auto& layout = self.GetLayout(); - const LO num_dof_holders = layout.GetNumOwnedDofHolder(); + const LO num_dof_holders = layout.GetNumDofHolder(); auto view = NumpyToHostComponentView(data, layout.GetNumComponents()); if (view.extent(0) != static_cast(num_dof_holders)) { @@ -135,7 +136,11 @@ void bind_create_field_module(py::module& m) py::arg("data"), "Set the DOF-holder values without checking their basis (the shape is " "still checked). A 1D array (num_dof_holders,) is accepted when the " - "field has one component."); + "field has one component.") + .def( + "synchronize", [](Field& self) { self.Synchronize(); }, + "Bring the copies of DOF holders other ranks own (e.g. Omega_h ghosts) " + "up to date with their owners' values"); py::class_, Field>( m, "Function", "A Field that knows its FunctionSpace") @@ -201,6 +206,9 @@ void bind_create_field_module(py::module& m) }, py::arg("name") = "", py::arg("variance") = values::Scalar, py::arg("system") = py::none(), py::arg("scaling") = py::none(), + // The function refers to the space's mesh (an Omega_h function keeps a + // tag on it), so the space, which keeps the mesh alive, outlives it. + py::keep_alive<0, 1>(), "Create a Function in this space. Vector and tensor values are in " "`system`'s basis (default: this space's coordinate system), with " "`scaling` required when that system admits both scalings.") diff --git a/src/pcms/pythonapi/uniform_grid_workflow_guide.md b/src/pcms/pythonapi/uniform_grid_workflow_guide.md index aca7bc2a..5befef7f 100644 --- a/src/pcms/pythonapi/uniform_grid_workflow_guide.md +++ b/src/pcms/pythonapi/uniform_grid_workflow_guide.md @@ -115,7 +115,7 @@ Create a mask where each uniform grid vertex is marked as 1 (inside mesh) or 0 ( ```python # Create binary mask indicating which vertices are inside the mesh mask_field = pcms.create_uniform_grid_binary_field(mesh, [4, 4]) -print(f"Created mask field with {mask_field.get_layout().get_num_owned_dof_holder()} vertices") +print(f"Created mask field with {mask_field.get_layout().get_num_dof_holder()} vertices") ``` **Accessing Values**: @@ -144,7 +144,7 @@ omega_h_field = omega_h_factory.create_function() # Initialize field with f(x,y) = x + 2*y coords = omega_h_field.get_layout().get_dof_holder_coordinates().get_values() -num_nodes = omega_h_field.get_layout().get_num_owned_dof_holder() +num_nodes = omega_h_field.get_layout().get_num_dof_holder() omega_h_data = np.zeros(num_nodes) for i in range(num_nodes): @@ -270,7 +270,7 @@ snake_case form of the C++ method (`GetDOFHolderData` is - `space.create_function(name="", variance=values.Scalar, system=None, scaling=None)` - Create a real-valued `Function`; vector and tensor values are in `system`'s basis (default: the space's) - `field.get_layout()` - the field's `FieldLayout` - `field.get_data().get_value_basis()` - the basis the field's values are stored in -- `field.get_layout().get_num_owned_dof_holder()` - Number of owned DOF holders (nodes/elements) +- `field.get_layout().get_num_dof_holder()` - Number of owned DOF holders (nodes/elements) - `field.get_layout().get_num_components()` - Number of field components per DOF holder - `field.get_layout().get_dof_holder_coordinates()` - DOF holder coordinates as a `CoordinateView` - `field.get_dof_holder_data()` - Field values as a `ValueView` whose array (`.get_values()`, shape `(num_dof_holders, num_components)`) is read-only and valid until the field is next modified (use `.copy()` to keep values) @@ -332,7 +332,7 @@ snake_case form of the C++ method (`GetDOFHolderData` is ### Grid Properties - `grid.get_num_cells()` - Total number of cells - `grid.divisions` - Cell divisions in each dimension -- `field.get_layout().get_num_owned_dof_holder()` - Total number of field DOF holders +- `field.get_layout().get_num_dof_holder()` - Total number of field DOF holders - `field.get_layout().get_dof_holder_coordinates().get_values()` - DOF holder coordinates --- diff --git a/src/pcms/transfer/copy.h b/src/pcms/transfer/copy.h index 73d63f32..528216a9 100644 --- a/src/pcms/transfer/copy.h +++ b/src/pcms/transfer/copy.h @@ -51,6 +51,7 @@ class Copy : public TransferOperator PCMS_FUNCTION_TIMER; detail::CheckCopyCompatible(source, target); target.SetDOFHolderData(source.GetDOFHolderData()); + target.Synchronize(); } }; diff --git a/src/pcms/transfer/interpolator.h b/src/pcms/transfer/interpolator.h index 959f5be3..f34420b0 100644 --- a/src/pcms/transfer/interpolator.h +++ b/src/pcms/transfer/interpolator.h @@ -67,6 +67,7 @@ class Interpolator : public TransferOperator evaluator_->Evaluate(source, ValueView( written, MakeRank2View(target_values_))); target.SetDOFHolderDataUnchecked(MakeConstRank2View(target_values_)); + target.Synchronize(); } private: diff --git a/src/pcms/transfer/omega_h_conservative_projection.cpp b/src/pcms/transfer/omega_h_conservative_projection.cpp index 45a0c1b6..a94145e1 100644 --- a/src/pcms/transfer/omega_h_conservative_projection.cpp +++ b/src/pcms/transfer/omega_h_conservative_projection.cpp @@ -58,8 +58,8 @@ OmegaHConservativeProjection::OmegaHConservativeProjection( solver_ = std::make_unique(mass_integrator, *rhs_integrator_); target_values_ = Kokkos::View( - "conservative_projection_target_values", - target_layout_->GetNumOwnedDofHolder(), target_layout_->GetNumComponents()); + "conservative_projection_target_values", target_layout_->GetNumDofHolder(), + target_layout_->GetNumComponents()); } // Defined here so that GalerkinProjectionSolver (forward-declared in the @@ -78,21 +78,25 @@ void OmegaHConservativeProjection::Apply(const Field& source, throw pcms_error( "OmegaHConservativeProjection::Apply: only scalar fields are supported"); } - const int num_dof_holders = target_layout_->GetNumOwnedDofHolder(); + const int num_dof_holders = target_layout_->GetNumDofHolder(); const int num_components = target_layout_->GetNumComponents(); const auto values = MakeRank2View(target_values_); const auto solution = solver_->Solve(*evaluator_, source); const auto global_to_local = target_layout_->GetGlobalToLocalPermutation(); + const bool has_ghosts = target_layout_->HasGhosts(); + const auto owned_indices = target_layout_->GetOwnedIndices(); Kokkos::parallel_for( "conservative_projection_scatter_solution", Kokkos::RangePolicy(0, num_dof_holders), - KOKKOS_LAMBDA(int i) { + KOKKOS_LAMBDA(int o) { + const LO i = has_ghosts ? owned_indices(o) : o; for (int c = 0; c < num_components; ++c) { - values(i, c) = solution[global_to_local(i) * num_components + c]; + values(o, c) = solution[global_to_local(i) * num_components + c]; } }); Kokkos::fence(); target.SetDOFHolderDataUnchecked(MakeConstRank2View(target_values_)); + target.Synchronize(); } } // namespace pcms diff --git a/src/pcms/transfer/omega_h_control_variate_projection.cpp b/src/pcms/transfer/omega_h_control_variate_projection.cpp index d8c4e1aa..123be95d 100644 --- a/src/pcms/transfer/omega_h_control_variate_projection.cpp +++ b/src/pcms/transfer/omega_h_control_variate_projection.cpp @@ -35,7 +35,7 @@ OmegaHControlVariateProjection::OmegaHControlVariateProjection( *rhs_integrator_); target_values_ = Kokkos::View( - "cv_target_values", target_layout_->GetNumOwnedDofHolder(), + "cv_target_values", target_layout_->GetNumDofHolder(), target_layout_->GetNumComponents()); const auto num_samples = sample_coords.GetValues().extent(0); f_samples_ = @@ -51,7 +51,7 @@ OmegaHControlVariateProjection::~OmegaHControlVariateProjection() = default; void OmegaHControlVariateProjection::Apply(const Field& source, Field& target) const { - const int num_dof_holders = target_layout_->GetNumOwnedDofHolder(); + const int num_dof_holders = target_layout_->GetNumDofHolder(); const int num_components = target_layout_->GetNumComponents(); // 1. Control variate: interpolate the source field onto the target space. @@ -79,10 +79,12 @@ void OmegaHControlVariateProjection::Apply(const Field& source, // 3. Project the residual: M * delta = r. const auto delta = solver_->Solve(MakeConstRank2View(residual)); - // 4. x = g + delta. delta is indexed by active PETSc row, while the field is - // indexed by local DOF holder. + // 4. x = g + delta. delta is indexed by active PETSc row, while the field + // values are indexed by owned DOF holder. const auto g_nodal = control_variate_.GetDOFHolderData().GetValues(); const auto global_to_local = target_layout_->GetGlobalToLocalPermutation(); + const bool has_ghosts = target_layout_->HasGhosts(); + const auto owned_indices = target_layout_->GetOwnedIndices(); if (static_cast(g_nodal.extent(0)) != num_dof_holders || static_cast(g_nodal.extent(1)) != num_components) { @@ -94,14 +96,16 @@ void OmegaHControlVariateProjection::Apply(const Field& source, Kokkos::parallel_for( "cv_add_correction", Kokkos::RangePolicy(0, num_dof_holders), - KOKKOS_LAMBDA(int i) { + KOKKOS_LAMBDA(int o) { + const LO i = has_ghosts ? owned_indices(o) : o; for (int c = 0; c < num_components; ++c) { - values(i, c) = - g_nodal(i, c) + delta[global_to_local(i) * num_components + c]; + values(o, c) = + g_nodal(o, c) + delta[global_to_local(i) * num_components + c]; } }); Kokkos::fence(); target.SetDOFHolderDataUnchecked(MakeConstRank2View(target_values_)); + target.Synchronize(); } } // namespace pcms diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index 974f3d0e..5e591f56 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -198,7 +198,7 @@ Data BuildDataImpl(const OmegaHLagrangeLayout& source_layout, d.coeffs = std::move(coeffs); d.ndof_per_elem = ndof; d.num_target_dofs = - static_cast(target_layout.GetNumOwnedDofHolder()); + static_cast(target_layout.GetNumDofHolderWithGhosts()); return d; } diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index c34dd529..0c897ba2 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -94,7 +94,7 @@ Mat BuildOmegaHMassMatrixImpl(Omega_h::Mesh& mesh, { const auto global_to_local = target_layout.GetGlobalToLocalPermutation(); const PetscInt num_dofs = - static_cast(target_layout.GetNumOwnedDofHolder()); + static_cast(target_layout.GetNumDofHolderWithGhosts()); const int nelems = mesh.nelems(); Mat mat = nullptr; diff --git a/src/pcms/utility/arrays.h b/src/pcms/utility/arrays.h index ca7f6490..4104603b 100644 --- a/src/pcms/utility/arrays.h +++ b/src/pcms/utility/arrays.h @@ -399,6 +399,29 @@ void CopyDeviceRank2ViewToDeviceView(Kokkos::View dest, }); } +// Inverse of CopyDeviceRank2ViewToDeviceView(Kokkos::View, Rank2View): +// scatter a flat node-major device View back into a [dof][comp] device +// Kokkos::View. T is deduced from the (non-alias) Kokkos::View +// destination parameter, since alias-template argument deduction through +// Rank2View's multi-level alias chain is not reliable. Uses 2D indexing on +// the destination so the result is correct regardless of its layout. +template +void CopyDeviceViewToRank2View(Kokkos::View dest, + Kokkos::View src) +{ + const size_t num_dof = dest.extent(0); + const size_t num_comp = dest.extent(1); + if (src.extent(0) != num_dof * num_comp) { + throw pcms_error("CopyDeviceViewToRank2View: size mismatch"); + } + Kokkos::parallel_for( + "CopyDeviceViewToRank2View", num_dof, KOKKOS_LAMBDA(LO i) { + for (size_t c = 0; c < num_comp; ++c) { + dest(i, c) = src(i * num_comp + c); + } + }); +} + // Host analogue: copy a [dof][comp] host Rank2View into a flat node-major // device View. The host Rank2View is layout_right (contiguous node-major), so // we alias its storage with an unmanaged view and deep_copy straight to device. @@ -502,5 +525,20 @@ void CopyRank1ViewToDevice(Rank1View dest, KOKKOS_LAMBDA(LO i) { dest(i) = src_tmp_d(i); }); } +// dest(rows(i), c) = src(i, c) for every row i of src, on the device. +template +void ScatterRows(Rank2View src, + Kokkos::View rows, + Rank2View dest) +{ + const size_t num_rows = src.extent(0); + const size_t num_comp = src.extent(1); + Kokkos::parallel_for( + "pcms::ScatterRows", + Kokkos::MDRangePolicy>( + {0, 0}, {num_rows, num_comp}), + KOKKOS_LAMBDA(size_t i, size_t c) { dest(rows(i), c) = src(i, c); }); +} + } // namespace pcms #endif // PCMS_COUPLING_ARRAYS_H diff --git a/src/pcms/utility/omega_h_array_utils.h b/src/pcms/utility/omega_h_array_utils.h index 9b76d97a..61bc3fa5 100644 --- a/src/pcms/utility/omega_h_array_utils.h +++ b/src/pcms/utility/omega_h_array_utils.h @@ -2,8 +2,11 @@ #define PCMS_UTILITY_OMEGA_H_ARRAY_UTILS_H #include +#include #include "pcms/utility/memory_spaces.h" +#include + namespace pcms { @@ -65,6 +68,46 @@ inline Kokkos::View ConvertCoordsTo2D( return coords_2d; } +// The value type Omega_h arrays and tags store for T. Omega_h supports I8, +// I32, I64 and Real, so a float is stored as Real. +template +struct OmegaHValue +{ + using type = T; +}; +template <> +struct OmegaHValue +{ + using type = Omega_h::Real; +}; +template +using OmegaHValue_t = typename OmegaHValue::type; + +// Synchronize a block of mesh entity data with Omega_h's internal arrays. +template +void SynchronizeOmegaHBlock(Omega_h::Mesh& mesh, int dim, int nc, + BlockView block) +{ + using Value = OmegaHValue_t; + if constexpr (std::is_same_v) { + auto synced = + mesh.sync_array(dim, Omega_h::Read(Omega_h::Write(block)), nc); + Kokkos::deep_copy(block, synced.view()); + } else { + // Omega_h has no array of T: synchronize a converted copy. + using exec = DeviceMemorySpace::execution_space; + const auto n = static_cast(block.size()); + Omega_h::Write arr(n, "field_data_block"); + Kokkos::parallel_for( + "pcms::SynchronizeOmegaHBlockIn", Kokkos::RangePolicy(0, n), + KOKKOS_LAMBDA(Omega_h::LO i) { arr[i] = static_cast(block(i)); }); + const auto synced = mesh.sync_array(dim, Omega_h::Read(arr), nc); + Kokkos::parallel_for( + "pcms::SynchronizeOmegaHBlockOut", Kokkos::RangePolicy(0, n), + KOKKOS_LAMBDA(Omega_h::LO i) { block(i) = static_cast(synced[i]); }); + } +} + } // namespace pcms #endif // PCMS_UTILITY_OMEGA_H_ARRAY_UTILS_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 97eee902..12ec4f11 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -34,6 +34,21 @@ if(PCMS_ENABLE_OMEGA_H) ${notRendezvous} ${cyclone1p}) endif() + add_exe(test_distributed_layout) + if(HOST_NPROC GREATER_EQUAL 2) + mpi_test(test_distributed_layout_2p 2 + $) + endif() + add_exe(test_distributed_field) + if(HOST_NPROC GREATER_EQUAL 2) + mpi_test(test_distributed_field_2p 2 + $) + endif() + add_exe(test_distributed_coupling) + if(HOST_NPROC GREATER_EQUAL 2) + mpi_test(test_distributed_coupling_2p 2 + $) + endif() add_exe(test_ohClassPtn_appRibPtn) set(cyclone2p ${PCMS_TEST_DATA_DIR}/cyclone/23elements/2p.osh/) if(HOST_NPROC GREATER_EQUAL 4) diff --git a/test/distributed_test_utils.h b/test/distributed_test_utils.h new file mode 100644 index 00000000..d1d92fa7 --- /dev/null +++ b/test/distributed_test_utils.h @@ -0,0 +1,147 @@ +#ifndef PCMS_TEST_DISTRIBUTED_TEST_UTILS_H +#define PCMS_TEST_DISTRIBUTED_TEST_UTILS_H + +// Shared helpers for the MPI tests of distributed fields. +// +// These tests run under mpirun without Catch2, so a failed check aborts every +// rank: a rank that just returned would leave the others blocked in the next +// collective. + +#include + +#include + +#include "pcms/configuration.h" +#include "pcms/field/data/omega_h.h" +#include "pcms/field/data/simple.h" +#ifdef PCMS_ENABLE_MESHFIELDS +#include "pcms/field/data/mesh_fields.h" +#endif +#include "pcms/field/field.h" +#include "pcms/field/value_view.hpp" + +#include +#include +#include +#include + +namespace distributed_test +{ + +using pcms::LO; +using pcms::Real; + +/// Abort every rank with `message` unless `condition` holds. +inline void require(bool condition, const char* message) +{ + if (!condition) { + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, "[rank %d] check failed: %s\n", rank, message); + MPI_Abort(MPI_COMM_WORLD, 1); + } +} + +/// Report `what` and abort every rank. For exceptions escaping `main`: a rank +/// that just returned would leave the others blocked in their next collective. +[[noreturn]] inline void AbortAllRanks(const char* what) noexcept +{ + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, "[rank %d] uncaught exception: %s\n", rank, what); + MPI_Abort(MPI_COMM_WORLD, 1); + std::abort(); // MPI_Abort does not return, but is not declared noreturn +} + +/// `values` tagged with `field`'s stored value basis: what a field's checked +/// setters and an evaluator's output expect. A Catch2-free counterpart of +/// `pcms::test::TagLike`. +template +pcms::ValueView TagLike( + const pcms::Field>& field, + pcms::Rank2View values) +{ + return {field.GetData().GetValueBasis(), values}; +} + +template +pcms::ValueView TagLike( + const pcms::Field& field, + Kokkos::View values) +{ + return {field.GetData().GetValueBasis(), pcms::MakeRank2View(values)}; +} + +/// The values of every DOF holder of `field`, owned and ghost, in mesh order, +/// copied to the host; ghost values are correct only after Synchronize. Field +/// itself only +/// exposes the owned values, so this reads the backend's storage: the Omega_h +/// tag of an OmegaHFieldData, or the MeshFields backend. +template +Kokkos::View ValuesWithGhostsHost( + const pcms::Field& field) +{ + if (const auto* omega_h_data = + dynamic_cast*>(&field.GetData())) { + const auto& layout = field.GetLayout(); + const auto nc = static_cast(layout.GetNumComponents()); + // The values of every DOF holder are in the field's Omega_h tag. + using TagValue = typename pcms::OmegaHFieldData::tag_value_type; + const auto& discretization = + dynamic_cast( + *layout.GetDiscretization()); + const auto tag = + Omega_h::HostRead(discretization.GetMesh().get_array( + layout.GetDOFHolderEntityDim(), omega_h_data->GetTagName())); + const size_t n = static_cast(tag.size()) / nc; + Kokkos::View host("values_with_ghosts_host", n, + nc); + for (size_t i = 0; i < n; ++i) { + for (size_t comp = 0; comp < nc; ++comp) { + host(i, comp) = static_cast(tag[i * nc + comp]); + } + } + return host; + } +#ifdef PCMS_ENABLE_MESHFIELDS + if (const auto* mesh_fields = + dynamic_cast*>(&field.GetData())) { + // The MeshFields backend stores every local entity; read it per entity + // dimension, in the layout's dimension-block order. + const auto& layout = + dynamic_cast(field.GetLayout()); + auto& mesh = layout.GetMesh(); + const auto nodes_per_dim = layout.GetNodesPerDim(); + const auto nc = static_cast(layout.GetNumComponents()); + Kokkos::View host( + "values_with_ghosts_host", + static_cast(layout.GetNumDofHolderWithGhosts()), nc); + size_t row = 0; + for (int dim = 0; dim <= mesh.dim(); ++dim) { + if (!nodes_per_dim[dim]) { + continue; + } + const size_t rows = static_cast(mesh.nents(dim)) * + static_cast(nodes_per_dim[dim]); + Kokkos::View flat("backend_flat", rows * nc); + mesh_fields->GetMeshFieldBackend()->GetData( + pcms::Rank1View(flat.data(), rows * nc), + nodes_per_dim[dim], nc, dim); + const auto flat_h = + Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), flat); + for (size_t r = 0; r < rows; ++r) { + for (size_t comp = 0; comp < nc; ++comp) { + host(row + r, comp) = flat_h(r * nc + comp); + } + } + row += rows; + } + return host; + } +#endif + throw std::runtime_error("ValuesWithGhostsHost: unsupported field data type"); +} + +} // namespace distributed_test + +#endif // PCMS_TEST_DISTRIBUTED_TEST_UTILS_H diff --git a/test/python/test_coordinate_transformations.py b/test/python/test_coordinate_transformations.py index 1db5c094..0b896ad1 100644 --- a/test/python/test_coordinate_transformations.py +++ b/test/python/test_coordinate_transformations.py @@ -206,21 +206,21 @@ def test_dof_holder_data_is_tagged(self, cylindrical_space): def test_set_checks_the_basis(self, cylindrical_space): b = _radial_field(cylindrical_space) - n = b.get_layout().get_num_owned_dof_holder() + n = b.get_layout().get_num_dof_holder() with pytest.raises(RuntimeError, match="does not match"): b.set_dof_holder_data( pcms.ValueView(_cartesian_vector(), np.zeros((n, 3)))) def test_set_checks_the_shape(self, cylindrical_space): b = _radial_field(cylindrical_space) - n = b.get_layout().get_num_owned_dof_holder() + n = b.get_layout().get_num_dof_holder() with pytest.raises(ValueError): b.set_dof_holder_data( pcms.ValueView(_cylindrical_vector(), np.zeros((n - 1, 3)))) def test_plain_arrays_only_for_scalars(self, cylindrical_space): b = _radial_field(cylindrical_space) - n = b.get_layout().get_num_owned_dof_holder() + n = b.get_layout().get_num_dof_holder() with pytest.raises(TypeError, match="need a ValueView"): b.set_dof_holder_data(np.zeros((n, 3))) b.set_dof_holder_data_unchecked(np.ones((n, 3))) diff --git a/test/python/test_field_copy.py b/test/python/test_field_copy.py index fa616aa2..02627309 100644 --- a/test/python/test_field_copy.py +++ b/test/python/test_field_copy.py @@ -35,7 +35,7 @@ def test_copy(self, world, dim, order, num_components): # Create original field and set data original = factory.create_function() layout = original.get_layout() - ndata = layout.get_num_owned_dof_holder() * layout.get_num_components() + ndata = layout.get_num_dof_holder() * layout.get_num_components() # Create sequential array of IDs ids = np.arange(ndata, dtype=np.float64).reshape( diff --git a/test/python/test_numpy_conversion.py b/test/python/test_numpy_conversion.py index 4d65a75d..ba81b191 100644 --- a/test/python/test_numpy_conversion.py +++ b/test/python/test_numpy_conversion.py @@ -34,7 +34,7 @@ def _coordinates(field): def _num_dof_holders(field): - return field.get_layout().get_num_owned_dof_holder() + return field.get_layout().get_num_dof_holder() def _linear_field(space): diff --git a/test/python/test_omega_h_field.py b/test/python/test_omega_h_field.py index 09d432aa..2fac8022 100644 --- a/test/python/test_omega_h_field.py +++ b/test/python/test_omega_h_field.py @@ -30,13 +30,13 @@ def test_field_methods(self, world, dim, order, num_components): field = factory.create_function() assert field.get_layout().get_num_components() == num_components - assert field.get_layout().get_num_owned_dof_holder() > 0 + assert field.get_layout().get_num_dof_holder() > 0 coords = field.get_layout().get_dof_holder_coordinates().get_values() - assert coords.shape[0] == field.get_layout().get_num_owned_dof_holder() + assert coords.shape[0] == field.get_layout().get_num_dof_holder() assert coords.shape[1] == dim - ndata = field.get_layout().get_num_owned_dof_holder() * num_components + ndata = field.get_layout().get_num_dof_holder() * num_components test_data = np.arange(ndata, dtype=np.float64).reshape( -1, num_components) expected_data = test_data.copy() @@ -62,9 +62,9 @@ def test_field_transfer(self, world, dim, order, num_components): target = target_space.create_function() coords = source.get_layout().get_dof_holder_coordinates().get_values() - num_owned = source.get_layout().get_num_owned_dof_holder() - source_data = np.zeros((num_owned, 1), dtype=np.float64) - for i in range(num_owned): + num_local = source.get_layout().get_num_dof_holder() + source_data = np.zeros((num_local, 1), dtype=np.float64) + for i in range(num_local): source_data[i, 0] = np.sum(coords[i, :dim]) source.set_dof_holder_data(source_data) @@ -86,14 +86,14 @@ def test_field_evaluation(self, world, dim, order, num_components): mesh_coords = mesh.coords() num_verts = mesh.nverts() - num_owned = field.get_layout().get_num_owned_dof_holder() - ndata = num_owned * num_components + num_local = field.get_layout().get_num_dof_holder() + ndata = num_local * num_components def test_func(x, y, z, component): return np.sin(5.0 * x * y) + float(component) test_data = np.zeros(ndata, dtype=np.float64) - for i in range(min(num_verts, num_owned)): + for i in range(min(num_verts, num_local)): x = mesh_coords[i * dim + 0] if dim >= 1 else 0.0 y = mesh_coords[i * dim + 1] if dim >= 2 else 0.0 z = mesh_coords[i * dim + 2] if dim >= 3 else 0.0 @@ -101,13 +101,13 @@ def test_func(x, y, z, component): idx = i * num_components + c test_data[idx] = test_func(x, y, z, c) - if order == 2 and num_owned > num_verts: - for i in range(num_verts, num_owned): + if order == 2 and num_local > num_verts: + for i in range(num_verts, num_local): for c in range(num_components): idx = i * num_components + c test_data[idx] = float(c) + 0.5 - field.set_dof_holder_data(test_data.reshape(num_owned, num_components)) + field.set_dof_holder_data(test_data.reshape(num_local, num_components)) if dim == 2: eval_coords = np.array([ diff --git a/test/python/test_uniform_grid_field.py b/test/python/test_uniform_grid_field.py index 286297be..b145a6ec 100644 --- a/test/python/test_uniform_grid_field.py +++ b/test/python/test_uniform_grid_field.py @@ -13,7 +13,7 @@ def _coordinates(field): def _num_dof_holders(field): - return field.get_layout().get_num_owned_dof_holder() + return field.get_layout().get_num_dof_holder() class TestUniformGridField: diff --git a/test/test_distributed_coupling.cpp b/test/test_distributed_coupling.cpp new file mode 100644 index 00000000..3d7530bb --- /dev/null +++ b/test/test_distributed_coupling.cpp @@ -0,0 +1,465 @@ +#include +#include +#include +#include +#include + +#include "pcms/coupler/field_exchange_planner.h" +#include "pcms/coupler/field_serializer.h" +#include "pcms/coupler/overlap_mask.h" +#include "pcms/field/data/simple.h" +#include "distributed_test_utils.h" +#include "pcms/field/coordinate_systems/cartesian.hpp" +#include "pcms/field/function_space/lagrange.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/types.h" + +#include +#include +#include +#include +#include + +using pcms::GO; +using pcms::LO; +using pcms::Real; + +namespace +{ +void require(bool cond, const char* msg) +{ + if (!cond) { + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, "[rank %d] distributed coupling check failed: %s\n", + rank, msg); + MPI_Abort(MPI_COMM_WORLD, 1); + } +} + +// Builds a one-block GID message advertising `gids`, i.e. the structure the +// exchange planner produces for a single destination. +std::vector MakeGidMessage(const std::vector& gids) +{ + require(!gids.empty(), "cannot build an empty GID message"); + std::vector msg( + static_cast(pcms::ent_offsets_len) + gids.size(), 0); + // Every DOF holder in these tests sits on a mesh vertex (entity dim 0), so + // the per-entity-dimension header counts equal the total GID count from slot + // one onwards. + for (int e = 1; e < pcms::ent_offsets_len; ++e) { + msg[static_cast(e)] = static_cast(gids.size()); + } + std::copy(gids.begin(), gids.end(), msg.begin() + pcms::ent_offsets_len); + return msg; +} + +// In-message layout for a one-block message (one sender group) received by +// `rank`. The planner reads only the entries describing this rank's segment of +// the receive buffer: its start and its end. +redev::InMessageLayout MakeInMessageLayout(int rank, int nproc, size_t msg_size) +{ + redev::InMessageLayout layout; + layout.srcRanks.assign(static_cast(nproc), 0); + layout.offset.assign(static_cast(nproc) + 1, 0); + layout.offset[static_cast(rank)] = 0; + layout.offset[static_cast(rank) + 1] = + static_cast(msg_size); + return layout; +} + +// --------------------------------------------------------------------------- +// Checks +// --------------------------------------------------------------------------- + +// Shared state for the distributed coupling checks: the (ghosted) mesh, its +// order-1 Lagrange space and the rank counts. Everything else the checks need +// (counts, GIDs, masks) is derived from the layout. +struct DistributedCouplingContext +{ + Omega_h::Mesh& mesh; + std::shared_ptr space; + int rank; + int nproc; + + const pcms::FieldLayout& Layout() const { return *space->GetLayout(); } + LO NumOwned() const { return Layout().GetNumDofHolder(); } + // Order-1 space: the DOF holders with ghosts are the mesh vertices. + LO NumLocal() const { return mesh.nverts(); } + int NumComponents() const { return Layout().GetNumComponents(); } +}; + +// Exercises the owned-only serialization contract used by distributed coupling: +// only owned (rank-exclusive) DOF holders are put on the wire, and received +// values are written back to the owned DOF holders of the field. +void check_serializer_round_trip(const DistributedCouplingContext& ctx) +{ + const pcms::FieldLayout& layout = ctx.Layout(); + const LO n_owned = ctx.NumOwned(); + const LO n_local = ctx.NumLocal(); + const LO n_comp = ctx.NumComponents(); + const auto owned_indices = layout.GetOwnedIndicesHost(); + require(n_owned > 0, "rank owns at least one DOF holder"); + require(n_owned >= 2, "rank owns at least two DOF holders"); + require(n_local > n_owned, "ghosted mesh has local > owned"); + + // Set the owned field data: value = index of the owned holder among the DOF + // holders with ghosts. + auto field = ctx.space->CreateFunction(); + std::vector owned_data(static_cast(n_owned)); + for (LO o = 0; o < n_owned; ++o) { + owned_data[static_cast(o)] = static_cast(owned_indices(o)); + } + field.SetDOFHolderDataHost(distributed_test::TagLike( + field, pcms::Rank2View( + owned_data.data(), n_owned, n_comp))); + + // Owned-indexed permutation: the last owned holder is "outside the overlap + // region" (perm = -1); the rest get compact buffer slots 0..n_participating. + const LO n_participating = n_owned - 1; + std::vector permutation(static_cast(n_owned), -1); + for (LO o = 0; o < n_participating; ++o) { + permutation[static_cast(o)] = o; + } + + // Serialize: only participating owned data is written to the buffer. + pcms::FieldSerializer serializer; + std::vector buffer(static_cast(n_participating) * n_comp); + const int sent = + serializer.Serialize(field, pcms::make_array_view(buffer), + pcms::make_const_array_view(permutation)); + require(sent == static_cast(n_participating) * n_comp, + "serialized size matches participating holders"); + + auto owned_field = + pcms::FlattenToRank1View(field.GetDOFHolderDataHost().GetValues()); + require(owned_field.size() == static_cast(n_owned), + "owned field size"); + for (LO o = 0; o < n_participating; ++o) { + require(buffer[static_cast(o)] == + owned_field[static_cast(o)], + "serialized value == owned value"); + } + + // Deserialize into a fresh field: participating owned values are restored; + // ghost and non-participating owned holders stay zero. + auto target = ctx.space->CreateFunction(); + serializer.Deserialize(target, pcms::make_const_array_view(buffer), + pcms::make_const_array_view(permutation)); + auto target_data = + pcms::FlattenToRank1View(target.GetDOFHolderDataHost().GetValues()); + require(target_data.size() == static_cast(n_owned), + "target owned size"); + + for (LO o = 0; o < n_participating; ++o) { + require(target_data[static_cast(o)] == + owned_data[static_cast(o)], + "participating owned value restored"); + } + require(target_data[static_cast(n_owned - 1)] == Real(0), + "non-participating owned holder is zero"); + + auto owned_mask = layout.GetOwnedHost(); + const auto target_local = distributed_test::ValuesWithGhostsHost(target); + require(target_local.extent(0) == static_cast(n_local), + "target local size"); + for (LO i = 0; i < n_local; ++i) { + if (!owned_mask[static_cast(i)]) { + require(target_local(i, 0) == Real(0), "ghost holder is zero"); + } + } +} + +// Exercises the distributed exchange end to end: an ExchangePlan packed into a +// GID message, GID matching on the receive side, and the payload round trip. +// The redev transport is replaced by MPI collectives (Alltoall/Alltoallv), and +// the exchange partners are the ranks that actually own this rank's ghost DOF +// holders, so the check holds for any partitioning and rank count. +void check_distributed_exchange(const DistributedCouplingContext& ctx) +{ + Omega_h::Mesh& mesh = ctx.mesh; + const pcms::FieldLayout& layout = ctx.Layout(); + const int rank = ctx.rank; + const int nproc = ctx.nproc; + const LO n_owned = ctx.NumOwned(); + const LO n_local = ctx.NumLocal(); + const LO n_comp = ctx.NumComponents(); + const auto owned_mask = layout.GetOwnedHost(); + auto owned_gids = layout.GetGidsHost(); + // Global IDs of the local (owned and ghost) vertices, from the mesh tag the + // layout was built from. + const auto local_gids = + Omega_h::HostRead(mesh.get_array(0, "global")); + auto owner_remotes = mesh.ask_owners(0); + auto owner_ranks = Omega_h::HostRead(owner_remotes.ranks); + + // The real exchange partners: this rank's ghost DOF holders, grouped by the + // rank that owns them. A ghost owned by rank n is a DOF this rank would have + // to obtain from n in a coupled exchange. + std::vector neighbor_set; + std::vector> ghosts_from(static_cast(nproc)); + { + std::set seen; + for (LO i = 0; i < n_local; ++i) { + if (!owned_mask[static_cast(i)]) { + const int owner = static_cast(owner_ranks[static_cast(i)]); + ghosts_from[static_cast(owner)].push_back( + static_cast(local_gids[i])); + seen.insert(owner); + } + } + neighbor_set.assign(seen.begin(), seen.end()); + } + require(!neighbor_set.empty(), + "the ghosted mesh should have at least one neighbour"); + + // Send-side plan: the planner is driven by a synthetic peer partition that + // routes every owned holder to one actual neighbour. The routing is synthetic + // (a single mesh has no peer partition to query), but it drives the same + // partition-based plan construction a coupled application uses. + const int destination = neighbor_set.front(); + std::vector ptn_ranks = {destination, rank}; + std::vector ptn_cuts = {0.0, 2.0}; // heap layout: cuts[1] is the cut + redev::RCBPtn peer_ptn(2, ptn_ranks, ptn_cuts); + redev::Partition exchange_partition{peer_ptn}; + + // Every owned holder participates; the unmapped (-1) path is exercised on the + // receive side, where the advertised list is a strict subset of this rank's + // owned holders. + pcms::OverlapMask overlap_mask(static_cast(n_owned)); + + pcms::GenericFieldExchangePlanner planner; + pcms::ExchangePlan send_plan = + planner.BuildExchangePlan(layout, exchange_partition, &overlap_mask); + + // Plan structure: a single destination (an actual neighbour) described by a + // valid CSR, and a permutation that compacts the owned holders in owned + // order. + require(send_plan.dest_ranks.size() == 1, + "send plan should route to a single destination"); + require(send_plan.dest_ranks[0] == destination, + "send plan destination should be an actual neighbour"); + require(send_plan.offsets.size() == 2, "send plan CSR size"); + require(send_plan.offsets[0] == 0, "send plan CSR start"); + require(send_plan.offsets[1] == static_cast(send_plan.msg_size), + "send plan CSR end != message size"); + require(send_plan.msg_size == static_cast(n_owned), + "every owned holder should participate"); + require(send_plan.permutation.size() == static_cast(n_owned), + "send plan permutation is owned-indexed"); + LO next_slot = 0; + for (LO o = 0; o < n_owned; ++o) { + require(send_plan.permutation[static_cast(o)] == next_slot, + "send plan permutation must compact owned holders in owned order"); + ++next_slot; + } + + // GID message: one block whose header counts the exchanged holders per mesh + // entity dimension and whose payload lists the owned GIDs. + const size_t gid_header = + send_plan.dest_ranks.size() * static_cast(pcms::ent_offsets_len); + std::vector gid_msg(send_plan.msg_size + gid_header, 0); + planner.FillGidMessage(layout, send_plan, + pcms::Rank1View( + gid_msg.data(), gid_msg.size())); + require(gid_msg[0] == 0, "GID message header must start at zero"); + for (int e = 1; e < pcms::ent_offsets_len; ++e) { + require(gid_msg[static_cast(e)] == + static_cast(send_plan.msg_size), + "GID message header must count the vertex holders"); + } + for (LO o = 0; o < n_owned; ++o) { + const auto slot = + static_cast(send_plan.permutation[static_cast(o)]); + require(gid_msg[gid_header + slot] == static_cast(owned_gids(o)), + "GID message payload != owned GIDs"); + } + + // Exchange each ghost GID with its owner. Two-phase all-to-all: first the + // per-neighbour counts, then the GID lists themselves. Afterwards ghost_recv + // holds, grouped by sender rank, every GID this rank owns that some neighbour + // admits holding -- i.e. the owned DOFs this rank must send values for. + std::vector send_counts(static_cast(nproc), 0); + std::vector recv_counts(static_cast(nproc), 0); + std::vector sdispls(static_cast(nproc), 0); + std::vector rdispls(static_cast(nproc), 0); + for (int n = 0; n < nproc; ++n) { + send_counts[static_cast(n)] = + static_cast(ghosts_from[static_cast(n)].size()); + } + MPI_Alltoall(send_counts.data(), 1, MPI_INT, recv_counts.data(), 1, MPI_INT, + MPI_COMM_WORLD); + std::vector ghost_send; + for (int n = 0; n < nproc; ++n) { + sdispls[static_cast(n)] = static_cast(ghost_send.size()); + ghost_send.insert(ghost_send.end(), + ghosts_from[static_cast(n)].begin(), + ghosts_from[static_cast(n)].end()); + } + LO total_recv = 0; + for (int n = 0; n < nproc; ++n) { + rdispls[static_cast(n)] = static_cast(total_recv); + total_recv += recv_counts[static_cast(n)]; + } + std::vector ghost_recv(static_cast(total_recv), 0); + MPI_Alltoallv(ghost_send.data(), send_counts.data(), sdispls.data(), + MPI_INT64_T, ghost_recv.data(), recv_counts.data(), + rdispls.data(), MPI_INT64_T, MPI_COMM_WORLD); + require(total_recv > 0, + "no neighbour admits holding any of this rank's owned DOFs"); + + // The received GID list is what the receive plan is built from: it gathers + // into a single block (one sender's message) every GID some neighbour admits + // holding, then matches those against the holders this rank owns. + const std::vector recv_gid_message = MakeGidMessage(ghost_recv); + pcms::ExchangePlan recv_plan = planner.BuildReceivePlan( + layout, + pcms::GlobalIDView(recv_gid_message.data(), + recv_gid_message.size()), + rank, nproc, MakeInMessageLayout(rank, nproc, recv_gid_message.size())); + + require(recv_plan.permutation.size() == static_cast(n_owned), + "receive plan permutation is owned-indexed"); + require(recv_plan.msg_size == static_cast(total_recv), + "receive plan payload length != advertised GID count"); + require(recv_plan.dest_ranks.size() == 1 && recv_plan.offsets.size() == 2, + "receive plan should describe a single sender"); + require(recv_plan.offsets[0] == 0 && + recv_plan.offsets[1] == static_cast(recv_plan.msg_size), + "receive plan CSR"); + // A rank's owned DOF can be a ghost on several neighbours (e.g. a corner + // vertex shared by four ranks), so the same GID may legitimately be + // advertised more than once. Matching only cares about the distinct set. + const std::set advertised_to_me(ghost_recv.begin(), ghost_recv.end()); + LO n_matched = 0; + for (LO o = 0; o < n_owned; ++o) { + const bool advertised = + advertised_to_me.count(static_cast(owned_gids(o))) > 0; + const LO slot = recv_plan.permutation[static_cast(o)]; + require((slot >= 0) == advertised, + "receive plan matched the wrong owned holders"); + if (advertised) { + ++n_matched; + require(slot < static_cast(recv_plan.msg_size), + "receive plan buffer index out of range"); + } + } + require(n_matched == static_cast(advertised_to_me.size()), + "receive plan must match every advertised GID against an owned " + "holder"); + require(n_matched < n_owned, + "every owned holder was advertised: the unmatched-permutation path " + "was not exercised"); + + // Local packing: Serialize this rank's owned values (the GIDs themselves) + // with the receive plan permutation and round-trip them through Deserialize + // to confirm each matched holder maps onto its slot (checked just below). The + // cross-rank broadcast after it sends the same values back to every + // advertiser via the transpose of the GID exchange. + auto exchange_field = ctx.space->CreateFunction(); + std::vector exchange_vals(static_cast(n_owned), Real(0)); + for (LO o = 0; o < n_owned; ++o) { + exchange_vals[static_cast(o)] = static_cast(owned_gids(o)); + } + exchange_field.SetDOFHolderDataHost(distributed_test::TagLike( + exchange_field, pcms::Rank2View( + exchange_vals.data(), n_owned, n_comp))); + + pcms::FieldSerializer exchange_serializer; + std::vector payload(recv_plan.msg_size, Real(0)); + const int packed = exchange_serializer.Serialize( + exchange_field, pcms::make_array_view(payload), + pcms::make_const_array_view(recv_plan.permutation)); + require(packed == static_cast(recv_plan.msg_size), + "packed payload size != receive plan message size"); + + // Local check of the packing: a matched holder's owned value lands at its + // plan slot, and holders the plan did not match stay untouched. + auto echo_field = ctx.space->CreateFunction(); + exchange_serializer.Deserialize( + echo_field, pcms::make_const_array_view(payload), + pcms::make_const_array_view(recv_plan.permutation)); + auto echo_data = + pcms::FlattenToRank1View(echo_field.GetDOFHolderDataHost().GetValues()); + for (LO o = 0; o < n_owned; ++o) { + if (recv_plan.permutation[static_cast(o)] >= 0) { + require(echo_data[static_cast(o)] == + static_cast(owned_gids(o)), + "packed value for a matched holder != its owned value"); + } else { + require(echo_data[static_cast(o)] == Real(0), + "unmatched holder must stay untouched by deserialization"); + } + } + const auto echo_local = distributed_test::ValuesWithGhostsHost(echo_field); + for (LO i = 0; i < n_local; ++i) { + if (!owned_mask[static_cast(i)]) { + require(echo_local(i, 0) == Real(0), + "deserialization must not write ghost holders"); + } + } + + // Cross-rank broadcast: the owner sends its value for every advertised GID + // back to each rank that advertised it, in that rank's own order. This is the + // same broadcast Omega_h's sync_array performs for ghost values -- it must + // serve every advertiser, because a corner DOF may be requested by several + // neighbours. + std::vector return_buf(ghost_recv.size(), Real(0)); + for (size_t k = 0; k < ghost_recv.size(); ++k) { + return_buf[k] = static_cast(ghost_recv[k]); + } + std::vector payload_recv(ghost_send.size(), Real(0)); + MPI_Alltoallv(return_buf.data(), recv_counts.data(), rdispls.data(), + MPI_DOUBLE, payload_recv.data(), send_counts.data(), + sdispls.data(), MPI_DOUBLE, MPI_COMM_WORLD); + for (size_t k = 0; k < ghost_send.size(); ++k) { + require(payload_recv[k] == static_cast(ghost_send[k]), + "values returned for ghost holders != the owning rank's values"); + } +} +} // namespace + +// Exercises the distributed coupling data path on a ghosted mesh: the +// owned-only serialization contract (check_serializer_round_trip) and the +// exchange plan, GID matching and payload round trip a coupler drives +// (check_distributed_exchange). +int main(int argc, char** argv) +{ + try { + MPI_Init(&argc, &argv); + int result = 0; + { + Omega_h::Library lib(&argc, &argv); + auto world = lib.world(); + const int rank = world->rank(); + const int nproc = world->size(); + require(nproc >= 2, "test requires at least 2 MPI ranks"); + + const int n = 4 * nproc; + Omega_h::Mesh mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, + 0.0, /*nx=*/n, /*ny=*/n, /*nz=*/0, + /*symmetric=*/false); + mesh.set_parting(OMEGA_H_GHOSTED, 1, false); + + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, /*order=*/1, /*num_components=*/1, + pcms::csys::Cartesian::Deferred(), "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + const DistributedCouplingContext ctx{mesh, space, rank, nproc}; + + check_serializer_round_trip(ctx); + + check_distributed_exchange(ctx); + + if (rank == 0) { + std::printf("distributed coupling test passed (nproc=%d)\n", nproc); + } + } + MPI_Finalize(); + return result; + } catch (const std::exception& e) { + distributed_test::AbortAllRanks(e.what()); + } catch (...) { + distributed_test::AbortAllRanks("unknown exception"); + } +} diff --git a/test/test_distributed_field.cpp b/test/test_distributed_field.cpp new file mode 100644 index 00000000..179a1947 --- /dev/null +++ b/test/test_distributed_field.cpp @@ -0,0 +1,252 @@ +#include + +#include +#include + +#include "distributed_test_utils.h" +#include "pcms/field/coordinate_systems/cartesian.hpp" +#include "pcms/field/function_space/lagrange.h" +#include "pcms/field/out_of_bounds_policy.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/types.h" + +#include +#include +#include +#include + +using pcms::LO; +using pcms::Real; + +namespace +{ +constexpr Real kFillValue = -2.0e30; + +// Sentinel used to pre-fill the evaluation buffer. It is deliberately different +// from kFillValue: the FILL policy writes kFillValue into every point that +// localization rejects, so if the buffer were pre-filled with the same value, a +// point the evaluator never wrote and a point it rejected would be +// indistinguishable. With a distinct sentinel the three outcomes separate: +// kUnresolved -> the evaluator never wrote this point (kernel/copy did not +// run) kFillValue -> localization rejected this point and the FILL policy +// wrote it anything else -> a real interpolated value +constexpr Real kUnresolved = -3.0e30; + +KOKKOS_INLINE_FUNCTION Real linear_f(Real x, Real y) +{ + return x + 2.0 * y; +} + +template +Kokkos::View copy_coordinates_to_host( + const CoordView& coordinates, LO n, int dim) +{ + using layout = std::conditional_t< + std::is_same_v, + Kokkos::LayoutLeft, Kokkos::LayoutRight>; + Kokkos::View + device_coordinates(coordinates.data_handle(), n, dim); + Kokkos::View host_coordinates("coordinates", n, + dim); + pcms::DeepCopyMismatchLayouts(host_coordinates, device_coordinates); + return host_coordinates; +} + +struct DeviceCoordinates +{ + Kokkos::View storage; + pcms::CoordinateView coordinate_view; +}; + +DeviceCoordinates make_device_coordinates(const std::vector& points) +{ + const LO n = static_cast(points.size() / 2); + Kokkos::View storage("coordinates", n, 2); + Kokkos::View host("coordinates", n, 2); + for (LO i = 0; i < n; ++i) { + host(i, 0) = points[static_cast(i) * 2]; + host(i, 1) = points[static_cast(i) * 2 + 1]; + } + pcms::DeepCopyMismatchLayouts(storage, host); + return {storage, + {pcms::csys::Cartesian::Deferred(), pcms::MakeRank2View(storage)}}; +} + +void require(bool condition, const char* message) +{ + if (!condition) { + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, "[rank %d] %s\n", rank, message); + MPI_Abort(MPI_COMM_WORLD, 1); + } +} + +void check_distributed_evaluation(Omega_h::Mesh& mesh, + const pcms::LagrangeFunctionSpace& space) +{ + const auto& layout = *space.GetLayout(); + // Order-1 layout: the DOF holders with ghosts are the mesh vertices. + const LO n_local = mesh.nverts(); + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + require(n_local > layout.GetNumDofHolder(), + "ghosted mesh did not create ghost DOF holders"); + + const LO n_owned = layout.GetNumDofHolder(); + // Coordinates of the mesh vertices, owned and ghost, one per DOF holder. + const auto mesh_coords = Omega_h::HostRead(mesh.coords()); + auto owned_coords = copy_coordinates_to_host( + layout.GetDOFHolderCoordinates().GetValues(), n_owned, mesh.dim()); + std::vector values(static_cast(n_owned)); + for (LO o = 0; o < n_owned; ++o) { + values[static_cast(o)] = + linear_f(owned_coords(o, 0), owned_coords(o, 1)); + } + + // Writing the owned values leaves the ghosts at their initial zero. + auto field = space.CreateFunction(); + field.SetDOFHolderDataHost(distributed_test::TagLike( + field, pcms::Rank2View(values.data(), + n_owned, 1))); + const auto owned = layout.GetOwnedHost(); + const auto before = distributed_test::ValuesWithGhostsHost(field); + for (LO i = 0; i < n_local; ++i) { + if (!owned[static_cast(i)]) { + require(before(i, 0) == 0.0, "owned write changed a ghost value"); + } + } + + field.Synchronize(); + const auto data = distributed_test::ValuesWithGhostsHost(field); + for (LO i = 0; i < n_local; ++i) { + const Real x = mesh_coords[i * mesh.dim()]; + const Real y = mesh_coords[i * mesh.dim() + 1]; + require(std::fabs(data(i, 0) - linear_f(x, y)) < 1.0e-12, + "ghost synchronization did not produce the affine field"); + } + + const LO n_elements = mesh.nents(mesh.dim()); + const int vertices_per_element = mesh.dim() + 1; + auto element_vertices = Omega_h::HostRead(mesh.ask_elem_verts()); + auto mesh_coordinates = Omega_h::HostRead(mesh.coords()); + std::vector points(static_cast(n_elements) * mesh.dim()); + for (LO element = 0; element < n_elements; ++element) { + for (int dim = 0; dim < mesh.dim(); ++dim) { + Real centroid = 0.0; + for (int vertex = 0; vertex < vertices_per_element; ++vertex) { + const LO vertex_id = + element_vertices[element * vertices_per_element + vertex]; + centroid += mesh_coordinates[vertex_id * mesh.dim() + dim]; + } + points[static_cast(element) * mesh.dim() + dim] = + centroid / vertices_per_element; + } + } + auto device_coords = make_device_coordinates(points); + auto evaluator = + space.CreatePointEvaluator(pcms::EvaluationRequest::FromCoordinates( + device_coords.coordinate_view, + pcms::OutOfBoundsPolicy{pcms::OutOfBoundsMode::FILL, kFillValue})); + + const LO n_points = n_elements; + Kokkos::View result("result", n_points, 1); + Kokkos::deep_copy(result, kUnresolved); + evaluator->Evaluate(field, distributed_test::TagLike(field, result)); + auto result_host = + Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), result); + + // A point still holding kUnresolved was never written by the evaluator (the + // evaluation kernels and/or the device-to-host copy produced nothing); a + // point holding kFillValue was written by the FILL policy, i.e. localization + // rejected it. They are counted separately because the causes differ. + LO n_unlocalized = 0; + LO n_untouched = 0; + LO first_unlocalized = -1; + LO first_untouched = -1; + for (LO i = 0; i < n_points; ++i) { + const Real value = result_host(i, 0); + if (value == kUnresolved) { + ++n_untouched; + if (first_untouched < 0) { + first_untouched = i; + } + continue; + } + if (value == kFillValue) { + ++n_unlocalized; + if (first_unlocalized < 0) { + first_unlocalized = i; + } + continue; + } + const size_t offset = static_cast(i) * 2; + const Real expected = linear_f(points[offset], points[offset + 1]); + require(std::fabs(value - expected) < 1.0e-12, + "distributed evaluation returned the wrong value"); + } + + if (n_untouched != 0) { + std::fprintf(stderr, + "[rank %d] %lld/%lld element centroid results were never " + "written by the evaluator (still at the pre-fill sentinel " + "%.17g, first index %lld): evaluation/copy failure, not a " + "localization failure\n", + rank, static_cast(n_untouched), + static_cast(n_points), + static_cast(kUnresolved), + static_cast(first_untouched)); + } + + if (n_unlocalized != 0) { + const size_t offset = static_cast(first_unlocalized) * 2; + std::fprintf(stderr, + "[rank %d] %lld/%lld element centroids were not localized " + "(first index %lld at point=(%.17g, %.17g))\n", + rank, static_cast(n_unlocalized), + static_cast(n_points), + static_cast(first_unlocalized), + static_cast(points[offset]), + static_cast(points[offset + 1])); + } + + // Distinct messages so the failure class is obvious from the abort line + // alone. + require(n_unlocalized == 0, "local element centroid was not localized"); + require(n_untouched == 0, + "element centroid evaluation result was never written"); + + if (rank == 0) + std::printf("distributed field test passed\n"); +} +} // namespace + +int main(int argc, char** argv) +{ + try { + MPI_Init(&argc, &argv); + { + Omega_h::Library library(&argc, &argv); + auto world = library.world(); + require(world->size() >= 2, "test requires at least 2 MPI ranks"); + + const int n = 4 * world->size(); + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, n, + n, 0, false); + mesh.set_parting(OMEGA_H_GHOSTED, 1, false); + auto space = pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, pcms::csys::Cartesian::Deferred(), "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH); + + check_distributed_evaluation(mesh, *space); + } + MPI_Finalize(); + return 0; + } catch (const std::exception& e) { + distributed_test::AbortAllRanks(e.what()); + } catch (...) { + distributed_test::AbortAllRanks("unknown exception"); + } +} diff --git a/test/test_distributed_layout.cpp b/test/test_distributed_layout.cpp new file mode 100644 index 00000000..f367e824 --- /dev/null +++ b/test/test_distributed_layout.cpp @@ -0,0 +1,271 @@ +// Checks the owned/ghost split of a FieldLayout on a ghosted, distributed +// Omega_h mesh: the DOF holders with ghosts are the owned ones plus their ghost +// copies, and the compact owned views +// (owned indices, owned GIDs, owned coordinates, owned field data) agree +// with the ones with ghosts. +#include + +#include +#include +#include +#include + +#include "distributed_test_utils.h" +#include "pcms/configuration.h" +#include "pcms/field/coordinate_systems/cartesian.hpp" +#include "pcms/field/function_space/lagrange.h" +#include "pcms/utility/arrays.h" +#include "pcms/utility/types.h" + +#include +#include +#include +#include +#include + +using pcms::GO; +using pcms::LO; +using pcms::Real; + +namespace +{ + +void require(bool condition, const std::string& message) +{ + if (!condition) { + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, "[rank %d] %s\n", rank, message.c_str()); + MPI_Abort(MPI_COMM_WORLD, 1); + } +} + +template +auto ToHost(const View& view) +{ + return Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), view); +} + +void CheckLayout(Omega_h::Mesh& mesh, const pcms::LagrangeFunctionSpace& space, + const std::string& label) +{ + const auto& layout = *space.GetLayout(); + // Order-1 layouts: the DOF holders with ghosts are the mesh vertices. + const LO n_local = mesh.nverts(); + const LO n_owned = layout.GetNumDofHolder(); + const int nc = layout.GetNumComponents(); + const int dim = layout.GetDimension(); + + require(layout.IsDistributed(), label + ": layout is not distributed"); + require(n_owned <= n_local, label + ": more owned than local holders"); + require(layout.Size() == n_owned * nc, label + ": Size mismatch"); + require(layout.HasGhosts() == (n_owned != n_local), + label + ": HasGhosts mismatch"); + + // Every holder is owned by exactly one rank, and a ghosted mesh has ghosts. + const long long n_owned_ll = n_owned; + const long long n_ghost_ll = n_local - n_owned; + long long owned_sum = 0; + long long ghost_sum = 0; + MPI_Allreduce(&n_owned_ll, &owned_sum, 1, MPI_LONG_LONG, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(&n_ghost_ll, &ghost_sum, 1, MPI_LONG_LONG, MPI_SUM, + MPI_COMM_WORLD); + require(owned_sum == static_cast(layout.GetNumGlobalDofHolder()), + label + ": owned counts do not sum to the global count"); + require(ghost_sum > 0, label + ": ghosted mesh produced no ghost holders"); + + // owned indices: one entry per owned holder, increasing, owned only. + const auto owned = layout.GetOwnedHost(); + const auto owned_indices = layout.GetOwnedIndicesHost(); + require(static_cast(owned_indices.size()) == n_owned, + label + ": owned->local map has the wrong size"); + LO n_owned_flags = 0; + for (LO i = 0; i < n_local; ++i) { + n_owned_flags += owned[i] ? 1 : 0; + } + require(n_owned_flags == n_owned, label + ": owned flags != owned count"); + for (LO o = 0; o < n_owned; ++o) { + const LO l = owned_indices(o); + require(l >= 0 && l < n_local && owned[l], + label + ": owned->local maps to a non-owned holder"); + require(o == 0 || owned_indices(o - 1) < l, + label + ": owned->local map is not increasing"); + } + const auto owned_to_local_device = ToHost(layout.GetOwnedIndices()); + require(static_cast(owned_to_local_device.size()) == n_owned, + label + ": device owned->local map has the wrong size"); + for (LO o = 0; o < n_owned; ++o) { + require(owned_to_local_device(o) == owned_indices(o), + label + ": device and host owned->local maps differ"); + } + + // Owned GIDs, coordinates and classification are the mesh's per-vertex + // arrays (ghosts included) at the mapped index. The layouts here are order-1, + // so DOF holder i is mesh vertex i. + const auto gids = + Omega_h::HostRead(mesh.get_array(0, "global")); + require(gids.size() == n_local, label + ": local count != mesh vertices"); + const auto owned_gids = layout.GetGidsHost(); + require(static_cast(owned_gids.size()) == n_owned, + label + ": owned GIDs have the wrong size"); + for (LO o = 0; o < n_owned; ++o) { + require(owned_gids[o] == gids[owned_indices(o)], + label + ": owned GID != local GID at the mapped index"); + } + + const auto coords = Omega_h::HostRead(mesh.coords()); + const auto owned_coords = layout.GetDOFHolderCoordinates().GetValues(); + require(static_cast(owned_coords.extent(0)) == n_owned, + label + ": owned coordinates have the wrong size"); + Kokkos::View owned_coords_copy( + "owned_coords", n_owned, dim); + pcms::CopyDeviceRank2ViewToRank2View(pcms::MakeRank2View(owned_coords_copy), + owned_coords); + const auto owned_coords_host = ToHost(owned_coords_copy); + for (LO o = 0; o < n_owned; ++o) { + for (int d = 0; d < dim; ++d) { + require(owned_coords_host(o, d) == coords[owned_indices(o) * dim + d], + label + ": owned coordinate != local coordinate"); + } + } + + // Owned entity offsets cover exactly the owned holders. + const auto owned_offsets = layout.GetEntOffsets(); + require(owned_offsets.back() == static_cast(n_owned), + label + ": owned entity offsets do not end at the owned count"); + + // Field data: write the owned values and read them back. The backend's owned + // entries are current right away and its ghosts keep their initial zero; + // after Synchronize, every DOF holder holds the value its owner wrote. + auto field = space.CreateFunction(); + std::vector values(static_cast(n_owned) * nc); + for (LO o = 0; o < n_owned; ++o) { + for (int c = 0; c < nc; ++c) { + values[static_cast(o) * nc + c] = + static_cast(owned_gids[o]) + 0.25 * c; + } + } + field.SetDOFHolderDataUncheckedHost( + pcms::Rank2View(values.data(), n_owned, + nc)); + + const auto owned_host = field.GetDOFHolderDataHost().GetValues(); + require(static_cast(owned_host.extent(0)) == n_owned && + static_cast(owned_host.extent(1)) == nc, + label + ": owned host data has the wrong shape"); + const auto owned_device = field.GetDOFHolderData().GetValues(); + require(static_cast(owned_device.extent(0)) == n_owned && + static_cast(owned_device.extent(1)) == nc, + label + ": owned device data has the wrong shape"); + Kokkos::View owned_device_copy( + "owned_device", n_owned, nc); + pcms::CopyDeviceRank2ViewToRank2View(pcms::MakeRank2View(owned_device_copy), + owned_device); + const auto owned_device_host = ToHost(owned_device_copy); + for (LO o = 0; o < n_owned; ++o) { + for (int c = 0; c < nc; ++c) { + const Real expected = static_cast(owned_gids[o]) + 0.25 * c; + require(owned_host(o, c) == expected, + label + ": owned host data != written value"); + require(owned_device_host(o, c) == expected, + label + ": owned device data != written value"); + } + } + + const auto before_sync = distributed_test::ValuesWithGhostsHost(field); + require(static_cast(before_sync.extent(0)) == n_local, + label + ": backend does not store every local holder"); + for (LO i = 0; i < n_local; ++i) { + for (int c = 0; c < nc; ++c) { + const Real expected = + owned[i] ? static_cast(gids[i]) + 0.25 * c : Real(0); + require(before_sync(i, c) == expected, + label + ": owned write missed an owned row or changed a ghost"); + } + } + + field.Synchronize(); + const auto local_values = distributed_test::ValuesWithGhostsHost(field); + for (LO i = 0; i < n_local; ++i) { + for (int c = 0; c < nc; ++c) { + const Real expected = static_cast(gids[i]) + 0.25 * c; + require(local_values(i, c) == expected, + label + ": a holder does not hold its owner's value after " + "Synchronize"); + } + } + + // Owned classification is the mesh classification at the mapped index. + const auto class_dims = + Omega_h::HostRead(mesh.get_array(0, "class_dim")); + const auto class_ids = Omega_h::HostRead( + mesh.get_array(0, "class_id")); + const auto owned_class_dims = + layout.GetDOFHolderClassificationDimensionsHost(); + const auto owned_class_ids = layout.GetDOFHolderClassificationIdsHost(); + require(static_cast(owned_class_dims.size()) == n_owned && + static_cast(owned_class_ids.size()) == n_owned, + label + ": owned classification has the wrong size"); + for (LO o = 0; o < n_owned; ++o) { + require(owned_class_dims[o] == class_dims[owned_indices(o)] && + owned_class_ids[o] == class_ids[owned_indices(o)], + label + ": owned classification != local classification"); + } +} + +// Report `what` and abort every rank: a rank that just returned would leave +// the others blocked in their next collective. +[[noreturn]] void AbortAllRanks(const char* what) noexcept +{ + int rank = -1; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + std::fprintf(stderr, "[rank %d] uncaught exception: %s\n", rank, what); + MPI_Abort(MPI_COMM_WORLD, 1); + std::abort(); // MPI_Abort does not return, but is not declared noreturn +} + +} // namespace + +int main(int argc, char** argv) +{ + try { + MPI_Init(&argc, &argv); + { + Omega_h::Library library(&argc, &argv); + auto world = library.world(); + require(world->size() >= 2, "test requires at least 2 MPI ranks"); + + const int n = 4 * world->size(); + auto mesh = Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, n, + n, 0, false); + mesh.set_parting(OMEGA_H_GHOSTED, 1, false); + + for (int nc : {1, 2}) { + const std::string suffix = " (" + std::to_string(nc) + " components)"; + CheckLayout(mesh, + *pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, nc, pcms::csys::Cartesian::Deferred(), "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH), + "Omega_h backend" + suffix); +#ifdef PCMS_ENABLE_MESHFIELDS + CheckLayout(mesh, + *pcms::LagrangeFunctionSpace::FromMesh( + mesh, 1, nc, pcms::csys::Cartesian::Deferred(), "global", + pcms::LagrangeFunctionSpace::Backend::MeshFields), + "MeshFields backend" + suffix); +#endif + } + + if (world->rank() == 0) { + std::printf("distributed layout test passed\n"); + } + } + MPI_Finalize(); + return 0; + } catch (const std::exception& e) { + AbortAllRanks(e.what()); + } catch (...) { + AbortAllRanks("unknown exception"); + } +} diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp index 7096c909..c01cb3be 100644 --- a/test/test_eqdsk.cpp +++ b/test/test_eqdsk.cpp @@ -188,7 +188,7 @@ TEST_CASE("EQDSKData with SplineFunctionSpace") auto [spline_space, psi_field] = pcms::MakeEQDSKField(eqdsk_data); auto layout = spline_space->GetLayout(); - const int layout_size = layout->OwnedSize(); + const int layout_size = layout->Size(); const int psizr_size = static_cast(eqdsk_data.PSIZR.extent(0)); REQUIRE(layout_size == psizr_size); diff --git a/test/test_evaluation_only_space.cpp b/test/test_evaluation_only_space.cpp index e76f2e1e..aff3778f 100644 --- a/test/test_evaluation_only_space.cpp +++ b/test/test_evaluation_only_space.cpp @@ -187,7 +187,7 @@ TEST_CASE("Functions of an evaluation-only space hold no values") { auto source = MockSpace::Create(1); auto f = source->CreateFunction("f", std::make_unique()); - REQUIRE(f.GetLayout().GetNumOwnedDofHolder() == 0); + REQUIRE(f.GetLayout().GetNumDofHolder() == 0); REQUIRE_THROWS_WITH(f.GetDOFHolderData(), ContainsSubstring("no DOFAccess")); REQUIRE_THROWS_WITH(source->CreateFunction("default"), ContainsSubstring("naming the quantity")); diff --git a/test/test_field_communication.cpp b/test/test_field_communication.cpp index 2a74e060..2922f753 100644 --- a/test/test_field_communication.cpp +++ b/test/test_field_communication.cpp @@ -124,7 +124,7 @@ void client1(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, mesh, order, 1, pcms::csys::Cartesian::Deferred()); auto layout = factory->GetLayout(); auto gids = layout->GetGidsHost(); - const auto n = layout->GetNumOwnedDofHolder(); + const auto n = layout->GetNumDofHolder(); Omega_h::HostWrite ids(n); PCMS_ALWAYS_ASSERT(n == gids.size()); Kokkos::parallel_for( @@ -156,7 +156,7 @@ void client2(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, mesh, order, 1, pcms::csys::Cartesian::Deferred()); auto layout = factory->GetLayout(); auto gids = layout->GetGidsHost(); - const auto n = layout->GetNumOwnedDofHolder(); + const auto n = layout->GetNumDofHolder(); auto field = factory->CreateFunction(); pcms::FieldLayoutCommunicator layout_comm(comm_name + "2", comm, rdv, channel, @@ -170,25 +170,20 @@ void client2(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, auto copied_array = pcms::FlattenToRank1View(field.GetDOFHolderDataHost().GetValues()); - auto owned = layout->GetOwnedHost(); PCMS_ALWAYS_ASSERT(copied_array.size() == gids.size()); - PCMS_ALWAYS_ASSERT(owned.size() == gids.size()); int expected = 0; int sum = 0; Kokkos::parallel_reduce( "reduce1", - Kokkos::RangePolicy(0, - owned.size()), - KOKKOS_LAMBDA(int i, int& local_sum) { local_sum += owned[i] != 0; }, - expected); + Kokkos::RangePolicy(0, gids.size()), + KOKKOS_LAMBDA(int /*i*/, int& local_sum) { local_sum += 1; }, expected); Kokkos::parallel_reduce( "reduce2", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(int i, int& local_sum) { - if (owned[i]) - local_sum += (int)copied_array[i] == gids[i]; + local_sum += (int)copied_array[i] == gids[i]; }, sum); @@ -216,7 +211,7 @@ void server(MPI_Comm comm, Omega_h::Mesh& mesh, std::string comm_name, auto factory = pcms::LagrangeFunctionSpace::FromMesh( mesh, order, 1, pcms::csys::Cartesian::Deferred()); auto layout = factory->GetLayout(); - const auto n = layout->GetNumOwnedDofHolder(); + const auto n = layout->GetNumDofHolder(); Omega_h::HostWrite ids(n); Kokkos::parallel_for( "id 0", Kokkos::RangePolicy(0, n), diff --git a/test/test_field_copy.cpp b/test/test_field_copy.cpp index 754371c6..29a8c2e0 100644 --- a/test/test_field_copy.cpp +++ b/test/test_field_copy.cpp @@ -24,7 +24,7 @@ void test_copy(Omega_h::CommPtr world, int dim, int order, int num_components) mesh, order, num_components, pcms::csys::Cartesian::Deferred(), "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); auto layout = factory->GetLayout(); - int ndata = layout->GetNumOwnedDofHolder() * num_components; + int ndata = layout->GetNumDofHolder() * num_components; Omega_h::HostWrite ids(ndata); Kokkos::parallel_for( Kokkos::RangePolicy(0, ndata), @@ -33,7 +33,7 @@ void test_copy(Omega_h::CommPtr world, int dim, int order, int num_components) auto original = factory->CreateFunction(); original.SetDOFHolderDataUncheckedHost( pcms::Rank2View( - ids.data(), layout->GetNumOwnedDofHolder(), num_components)); + ids.data(), layout->GetNumDofHolder(), num_components)); auto copied = factory->CreateFunction(); pcms::Copy copy(*factory, *factory); diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp index 85a05a96..d3ded87f 100644 --- a/test/test_field_evaluation.cpp +++ b/test/test_field_evaluation.cpp @@ -72,7 +72,8 @@ TEST_CASE("evaluate quadratic 2d meshfields_field") // per entity dimension in the same block order the layout stores: // dim 0 (vertex) holders -> the vertex's own coordinate, // dim 1 (edge) holders -> the midpoint of the edge's two endpoint vertices. -// Returns a flat node-major vector matching GetDOFHolderCoordinates() layout. +// Returns a flat node-major vector matching GetDOFHolderCoordinates() +// layout. static std::vector ComputeManualHolderCoords( Omega_h::Mesh& mesh, const std::array& nodes_per_dim) { diff --git a/test/test_field_serializer.cpp b/test/test_field_serializer.cpp index 26756492..9fbe83fa 100644 --- a/test/test_field_serializer.cpp +++ b/test/test_field_serializer.cpp @@ -211,7 +211,7 @@ TEST_CASE("FieldSerializer round-trips a multi-component Omega_h field") pcms::LagrangeFunctionSpace::Backend::OmegaH); auto field = space->CreateFunction(); - const int num_dof = field.GetLayout().GetNumOwnedDofHolder(); + const int num_dof = field.GetLayout().GetNumDofHolder(); std::vector data(static_cast(num_dof) * num_components); for (int i = 0; i < num_dof; ++i) { for (int c = 0; c < num_components; ++c) { diff --git a/test/test_mfem_adapter.cpp b/test/test_mfem_adapter.cpp index ba5b869a..87e43015 100644 --- a/test/test_mfem_adapter.cpp +++ b/test/test_mfem_adapter.cpp @@ -35,7 +35,7 @@ TEST_CASE("MFEM vertex-scalar field adapter", "[mfem]") REQUIRE(layout.GetNumComponents() == 1); REQUIRE(layout.GetDimension() == 2); - REQUIRE(layout.GetNumOwnedDofHolder() == pmesh.GetNV()); + REQUIRE(layout.GetNumDofHolder() == pmesh.GetNV()); REQUIRE(layout.IsDistributed()); auto gids = layout.GetGidsHost(); @@ -114,7 +114,7 @@ TEST_CASE("MFEM vertex-scalar field adapter", "[mfem]") auto field = factory.CreateField("source"); const auto& layout = field.GetLayout(); - const auto n = static_cast(layout.GetNumOwnedDofHolder()); + const auto n = static_cast(layout.GetNumDofHolder()); // Identity permutation: buffer index equals DOF-holder index. Kokkos::View perm("perm", n); diff --git a/test/test_mfem_coupling.cpp b/test/test_mfem_coupling.cpp index 495495af..8e80cfa1 100644 --- a/test/test_mfem_coupling.cpp +++ b/test/test_mfem_coupling.cpp @@ -85,13 +85,13 @@ int RunClient(MPI_Comm comm) app->SetLayoutOverlapMask( LayoutName, std::make_unique( - static_cast(layout->GetNumOwnedDofHolder()), overlap_view)); + static_cast(layout->GetNumDofHolder()), overlap_view)); auto handle = app->AddField(factory.CreateField(FieldName)); // Seed the field so each vertex holds (gid + 1). auto gids = layout->GetGidsHost(); - std::vector values(static_cast(layout->GetNumOwnedDofHolder())); + std::vector values(static_cast(layout->GetNumDofHolder())); for (size_t v = 0; v < values.size(); ++v) { values[v] = ExpectedValue(gids[v]); } @@ -119,10 +119,9 @@ int RunServer(MPI_Comm comm) auto layout = factory.GetLayout(); auto handle = app->AddField(factory.CreateField(FieldName)); - FillField( - handle.GetField(), - std::vector(static_cast(layout->GetNumOwnedDofHolder()), - OutsideOverlapSentinel)); + FillField(handle.GetField(), + std::vector(static_cast(layout->GetNumDofHolder()), + OutsideOverlapSentinel)); app->ReceivePhase([&]() { handle.Receive(); }); diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp index 65871ac5..2ae8df4d 100644 --- a/test/test_omega_h_lagrange_field.cpp +++ b/test/test_omega_h_lagrange_field.cpp @@ -36,7 +36,7 @@ TEST_CASE("OmegaHLagrangeLayout order-1 properties") REQUIRE(layout.GetOrder() == 1); REQUIRE(layout.GetNumComponents() == 2); - REQUIRE(layout.GetNumOwnedDofHolder() == mesh.nents(0)); + REQUIRE(layout.GetNumDofHolder() == mesh.nents(0)); REQUIRE(layout.GetNumGlobalDofHolder() == mesh.nglobal_ents(0)); REQUIRE(layout.IsDistributed()); // always true for Omega_h mesh layouts @@ -70,7 +70,7 @@ TEST_CASE("OmegaHLagrangeLayout order-0 properties") REQUIRE(layout.GetOrder() == 0); REQUIRE(layout.GetNumComponents() == 1); - REQUIRE(layout.GetNumOwnedDofHolder() == mesh.nelems()); + REQUIRE(layout.GetNumDofHolder() == mesh.nelems()); REQUIRE(layout.GetNumGlobalDofHolder() == mesh.nglobal_ents(mesh.dim())); // DOF holder coordinates should match element centroids @@ -130,7 +130,7 @@ TEST_CASE("OmegaHLagrangeField order-1: set/get DOF data round-trip") mesh, 1, 1, pcms::csys::Cartesian::Deferred()); auto field = factory->CreateFunction(); - int n = factory->GetLayout()->GetNumOwnedDofHolder(); + int n = factory->GetLayout()->GetNumDofHolder(); std::vector data(n); for (int i = 0; i < n; ++i) data[i] = static_cast(i); @@ -154,7 +154,7 @@ TEST_CASE("Field checked writes validate the view shape, not just total size") mesh, 1, 1, pcms::csys::Cartesian::Deferred()); auto field = factory->CreateFunction(); - const int n = factory->GetLayout()->GetNumOwnedDofHolder(); + const int n = factory->GetLayout()->GetNumDofHolder(); std::vector data(static_cast(n), 1.0); const auto& basis = field.GetData().GetValueBasis(); @@ -235,7 +235,7 @@ TEST_CASE("OmegaHLagrangeField order-0: set/get DOF data round-trip") mesh, 0, 1, pcms::csys::Cartesian::Deferred()); auto field = factory->CreateFunction(); - int n = factory->GetLayout()->GetNumOwnedDofHolder(); + int n = factory->GetLayout()->GetNumDofHolder(); std::vector data(n, 3.14); pcms::Rank2View view(data.data(), n, 1); field.SetDOFHolderDataHost( diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp index 02f2e20b..5adfe384 100644 --- a/test/test_point_evaluator.cpp +++ b/test/test_point_evaluator.cpp @@ -405,7 +405,7 @@ TEST_CASE("SimpleFieldData: set and get DOF holder data round-trip") auto field_data = factory->CreateFunction(); auto& layout = *factory->GetLayout(); - int n = layout.GetNumOwnedDofHolder(); + int n = layout.GetNumDofHolder(); REQUIRE(n > 0); // Write sequential values @@ -547,7 +547,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 multi-component (2) evaluation") auto field = factory->CreateFunction(); auto layout = factory->GetLayout(); REQUIRE(layout->GetNumComponents() == 2); - int num_dof = layout->GetNumOwnedDofHolder(); + int num_dof = layout->GetNumDofHolder(); // Component 0: f(x,y)=x+y, Component 1: f(x,y)=2x-y auto dof_coords_mdspan = layout->GetDOFHolderCoordinates().GetValues(); @@ -602,7 +602,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 multi-component (3) evaluation") auto field = factory->CreateFunction(); auto layout = factory->GetLayout(); REQUIRE(layout->GetNumComponents() == 3); - int num_dof = layout->GetNumOwnedDofHolder(); + int num_dof = layout->GetNumDofHolder(); // Component 0: x, Component 1: y, Component 2: x*y auto dof_coords_mdspan = layout->GetDOFHolderCoordinates().GetValues(); diff --git a/test/test_polynomial_reconstruction_function_space.cpp b/test/test_polynomial_reconstruction_function_space.cpp index c6286529..bb0148f5 100644 --- a/test/test_polynomial_reconstruction_function_space.cpp +++ b/test/test_polynomial_reconstruction_function_space.cpp @@ -42,7 +42,7 @@ TEST_CASE( auto layout = factory->GetLayout(); REQUIRE(layout->GetNumComponents() == 1); - REQUIRE(layout->GetNumOwnedDofHolder() == 4); + REQUIRE(layout->GetNumDofHolder() == 4); REQUIRE(layout->GetNumGlobalDofHolder() == 4); REQUIRE_FALSE(layout->IsDistributed()); @@ -114,7 +114,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace field keeps layout alive " auto point_cloud_layout = dynamic_cast(&field.GetLayout()); REQUIRE(point_cloud_layout != nullptr); - REQUIRE(point_cloud_layout->GetNumOwnedDofHolder() == 4); + REQUIRE(point_cloud_layout->GetNumDofHolder() == 4); std::vector data{9.0, 10.0, 11.0, 12.0}; Rank2View data_view( @@ -144,8 +144,8 @@ TEST_CASE("Different layouts on the same mesh report SameEntities") REQUIRE(lagrange_disc != nullptr); REQUIRE(nodal_disc->SameEntities(*lagrange_disc)); - REQUIRE(nodal->GetLayout()->GetNumOwnedDofHolder() == mesh.nfaces()); - REQUIRE(lagrange->GetLayout()->GetNumOwnedDofHolder() == mesh.nverts()); + REQUIRE(nodal->GetLayout()->GetNumDofHolder() == mesh.nfaces()); + REQUIRE(lagrange->GetLayout()->GetNumDofHolder() == mesh.nverts()); REQUIRE(nodal_disc->GetNumEntities(pcms::Face) == mesh.nfaces()); REQUIRE(lagrange_disc->GetNumEntities(pcms::Vertex) == mesh.nverts()); diff --git a/test/test_proxy_coupling.cpp b/test/test_proxy_coupling.cpp index 0e8c04d4..170bfee0 100644 --- a/test/test_proxy_coupling.cpp +++ b/test/test_proxy_coupling.cpp @@ -27,7 +27,7 @@ void initializeFieldWithGids(const pcms::FieldLayout& layout, pcms::Real multiplier = 1.0) { auto gids = layout.GetGidsHost(); - const auto n = layout.GetNumOwnedDofHolder(); + const auto n = layout.GetNumDofHolder(); Omega_h::HostWrite ids(n); PCMS_ALWAYS_ASSERT(n == gids.size()); @@ -51,26 +51,21 @@ bool validateField(const pcms::FieldLayout& layout, auto gids = layout.GetGidsHost(); auto copied_array = pcms::FlattenToRank1View(field->GetDOFHolderDataHost().GetValues()); - auto owned = layout.GetOwnedHost(); - const auto n = layout.GetNumOwnedDofHolder(); + const auto n = layout.GetNumDofHolder(); PCMS_ALWAYS_ASSERT(copied_array.size() == gids.size()); - PCMS_ALWAYS_ASSERT(owned.size() == gids.size()); int expected = 0; int sum = 0; Kokkos::parallel_reduce( "count_owned", - Kokkos::RangePolicy(0, - owned.size()), - KOKKOS_LAMBDA(int i, int& local_sum) { local_sum += owned[i] != 0; }, - expected); + Kokkos::RangePolicy(0, gids.size()), + KOKKOS_LAMBDA(int /*i*/, int& local_sum) { local_sum += 1; }, expected); Kokkos::parallel_reduce( "count_match", Kokkos::RangePolicy(0, n), KOKKOS_LAMBDA(int i, int& local_sum) { - if (owned[i]) - local_sum += (int)copied_array[i] == (int)(gids[i] * multiplier); + local_sum += (int)copied_array[i] == (int)(gids[i] * multiplier); }, sum); diff --git a/test/test_proxy_coupling_xgc_server.cpp b/test/test_proxy_coupling_xgc_server.cpp index b7d06754..0dcad6a7 100644 --- a/test/test_proxy_coupling_xgc_server.cpp +++ b/test/test_proxy_coupling_xgc_server.cpp @@ -3,7 +3,7 @@ #include #include #include -#include +#include #include "test_support.h" #include "pcms/field/function_space/xgc.h" @@ -174,7 +174,7 @@ void omegah_coupler(MPI_Comm comm, Omega_h::Mesh& mesh, mesh, 1, 1, pcms::csys::Cartesian::Deferred(), numbering, pcms::LagrangeFunctionSpace::Backend::OmegaH, ss.str()); auto field = factory->CreateFunction( - ss.str(), std::make_unique>( + ss.str(), std::make_unique>( factory->GetLayout(), pcms::ValueBasis{})); fields.push_back(application->AddField(std::move(field))); } diff --git a/test/test_proxy_coupling_xgc_server_overlap.cpp b/test/test_proxy_coupling_xgc_server_overlap.cpp index 4fd3eba1..7e798cec 100644 --- a/test/test_proxy_coupling_xgc_server_overlap.cpp +++ b/test/test_proxy_coupling_xgc_server_overlap.cpp @@ -3,7 +3,7 @@ #include #include #include -#include +#include #include "test_support.h" #include "pcms/field/function_space/xgc.h" @@ -174,7 +174,7 @@ void omegah_coupler_with_overlap(MPI_Comm comm, Omega_h::Mesh& mesh, pcms::LagrangeFunctionSpace::Backend::OmegaH, ss.str()); auto field = factory->CreateFunction( - ss.str(), std::make_unique>( + ss.str(), std::make_unique>( factory->GetLayout(), pcms::ValueBasis{})); fields.push_back(application->AddField(std::move(field))); diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index e4e44ec6..b5bd2238 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -105,7 +105,7 @@ TEST_CASE("UniformGrid field creation") grid, 1, pcms::csys::Cartesian::Deferred()); REQUIRE(layout->GetNumComponents() == 1); - REQUIRE(layout->GetNumOwnedDofHolder() == 36); // (5+1)x(5+1) = 36 vertices + REQUIRE(layout->GetNumDofHolder() == 36); // (5+1)x(5+1) = 36 vertices REQUIRE(layout->GetNumGlobalDofHolder() == 36); REQUIRE_FALSE(layout->IsDistributed()); @@ -113,7 +113,7 @@ TEST_CASE("UniformGrid field creation") grid, 1, pcms::csys::Cartesian::Deferred()); auto field = field_space->CreateFunction(); REQUIRE(field.GetDOFHolderDataHost().size() == - static_cast(layout->OwnedSize())); + static_cast(layout->Size())); } TEST_CASE("UniformGrid order-0 field creation and evaluation") @@ -131,7 +131,7 @@ TEST_CASE("UniformGrid order-0 field creation and evaluation") pcms::UniformGridEvaluatorFactory<2> eval_factory(layout); REQUIRE(layout->GetOrder() == 0); - REQUIRE(layout->GetNumOwnedDofHolder() == 4); + REQUIRE(layout->GetNumDofHolder() == 4); auto coords_device = layout->GetDOFHolderCoordinates().GetValues(); auto coords = pcms::test::CopyCoordinatesToHost(coords_device); @@ -316,7 +316,7 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") auto transferred_data = pcms::FlattenToRank1View(ug_field.GetDOFHolderDataHost().GetValues()); auto ug_coords = ug_factory->GetLayout()->GetDOFHolderCoordinates(); - int num_ug_nodes = ug_factory->GetLayout()->GetNumOwnedDofHolder(); + int num_ug_nodes = ug_factory->GetLayout()->GetNumDofHolder(); // set up_coords to host auto ug_coords_host = diff --git a/test/testing.cmake b/test/testing.cmake index d26a4a4e..a1e7e94f 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -102,12 +102,12 @@ function(mpi_test TESTNAME PROCS EXE) if(NOT ${MPIEXEC_PREFLAGS} MATCHES "none") list(APPEND TEST_COMMAND ${MPIEXEC_PREFLAGS}) endif() - list(APPEND TEST_COMMAND ARGS ${MPIEXEC_NUMPROC_FLAG} ${PROCS}) + list(APPEND TEST_COMMAND ${MPIEXEC_NUMPROC_FLAG} ${PROCS}) if(NOT ${VALGRIND_EXECUTABLE} MATCHES "none") list(APPEND TEST_COMMAND ${VALGRIND_EXECUTABLE} ${VALGRIND_ARGS}) endif() list(APPEND TEST_COMMAND ${EXE} ${ARGN}) - add_test(NAME ${TESTNAME} COMMAND ${TEST_COMMAND }) + add_test(NAME ${TESTNAME} COMMAND ${TEST_COMMAND}) set_mpi_test_environment(${TESTNAME}) isolate_mpi_test(${TESTNAME} ${PROCS}) endfunction(mpi_test)