diff --git a/include/compute/receivers/receivers.hpp b/include/compute/receivers/receivers.hpp index ca981288..9fcb7d83 100644 --- a/include/compute/receivers/receivers.hpp +++ b/include/compute/receivers/receivers.hpp @@ -329,9 +329,10 @@ struct receivers : public impl::StationIterator, * @return Kokkos::View View of the * elements indices associated with the receivers */ - Kokkos::View - get_elements_on_device(const specfem::element::medium_tag medium, - const specfem::element::property_tag property) const; + std::tuple, + Kokkos::View > + get_indices_on_device(const specfem::element::medium_tag medium, + const specfem::element::property_tag property) const; /** * @brief Get the spectral element indices in which the receivers are located @@ -344,9 +345,10 @@ struct receivers : public impl::StationIterator, * @return Kokkos::View View of the * elements indices associated with the receivers */ - Kokkos::View - get_elements_on_host(const specfem::element::medium_tag medium, - const specfem::element::property_tag property) const; + std::tuple, + Kokkos::View > + get_indices_on_host(const specfem::element::medium_tag medium, + const specfem::element::property_tag property) const; /** * @brief Get the seismogram types @@ -364,8 +366,6 @@ struct receivers : public impl::StationIterator, IndexViewType::HostMirror h_elements; ///< Host view to store the ///< elements associated with the ///< receivers - IndexViewType receiver_domain_index_mapping; - IndexViewType::HostMirror h_receiver_domain_index_mapping; LagrangeInterpolantType lagrange_interpolant; ///< Lagrange interpolant for ///< every receiver LagrangeInterpolantType::HostMirror @@ -380,6 +380,12 @@ struct receivers : public impl::StationIterator, GET_NAME(PROPERTY_TAG)); \ IndexViewType::HostMirror CREATE_VARIABLE_NAME( \ h_elements, GET_NAME(DIMENSION_TAG), GET_NAME(MEDIUM_TAG), \ + GET_NAME(PROPERTY_TAG)); \ + IndexViewType CREATE_VARIABLE_NAME( \ + receiver_indices, GET_NAME(DIMENSION_TAG), GET_NAME(MEDIUM_TAG), \ + GET_NAME(PROPERTY_TAG)); \ + IndexViewType::HostMirror CREATE_VARIABLE_NAME( \ + h_receiver_indices, GET_NAME(DIMENSION_TAG), GET_NAME(MEDIUM_TAG), \ GET_NAME(PROPERTY_TAG)); CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS( @@ -439,16 +445,9 @@ load_on_device(const MemberType &team_member, const IteratorType &iterator, Kokkos::abort(message.c_str()); } - if (receivers.receiver_domain_index_mapping(index.ispec) == -1) { - std::string message = "Invalid element detected in kernel at " + - std::string(__FILE__) + ":" + - std::to_string(__LINE__); - Kokkos::abort(message.c_str()); - } - #endif - const int irec = receivers.receiver_domain_index_mapping(index.ispec); + const int irec = iterator_index.imap; lagrange_interpolant(iterator_index.ielement, index.iz, index.ix, 0) = receivers.lagrange_interpolant(irec, index.iz, index.ix, 0); @@ -500,16 +499,9 @@ store_on_device(const MemberType &team_member, const IteratorType &iterator, Kokkos::abort(message.c_str()); } - if (receivers.receiver_domain_index_mapping(index.ispec) == -1) { - std::string message = "Invalid element detected in kernel at " + - std::string(__FILE__) + - std::to_string(__LINE__); - Kokkos::abort(message.c_str()); - } - #endif - const int irec = receivers.receiver_domain_index_mapping(index.ispec); + const int irec = iterator_index.imap; receivers.seismogram_components(isig_step, iseis, irec, 0) = seismogram_components(iterator_index.ielement, 0); diff --git a/include/kokkos_kernels/impl/compute_seismogram.tpp b/include/kokkos_kernels/impl/compute_seismogram.tpp index d0b2405c..57a5cbea 100644 --- a/include/kokkos_kernels/impl/compute_seismogram.tpp +++ b/include/kokkos_kernels/impl/compute_seismogram.tpp @@ -27,12 +27,12 @@ void specfem::kokkos_kernels::impl::compute_seismograms( constexpr int ngll = NGLL; constexpr auto dimension = DimensionType; - const auto elements = - assembly.receivers.get_elements_on_device(MediumTag, PropertyTag); + const auto [elements, receiver_indices] = + assembly.receivers.get_indices_on_device(MediumTag, PropertyTag); - const int nelements = elements.extent(0); + const int nreceivers = receiver_indices.extent(0); - if (nelements == 0) + if (nreceivers == 0) return; auto &receivers = assembly.receivers; @@ -61,7 +61,7 @@ void specfem::kokkos_kernels::impl::compute_seismograms( lane_size, no_simd, Kokkos::DefaultExecutionSpace>; - using ChunkPolicy = specfem::policy::element_chunk; + using ChunkPolicy = specfem::policy::mapped_element_chunk; using ChunkElementFieldType = specfem::chunk_element::field< ParallelConfig::chunk_size, ngll, DimensionType, MediumTag, specfem::kokkos::DevScratchSpace, Kokkos::MemoryTraits, @@ -84,7 +84,7 @@ void specfem::kokkos_kernels::impl::compute_seismograms( receivers.set_seismogram_step(isig_step); - ChunkPolicy policy(elements, ngll, ngll); + ChunkPolicy policy(elements, receiver_indices, ngll, ngll); for (int iseis = 0; iseis < nseismograms; ++iseis) { @@ -112,12 +112,12 @@ void specfem::kokkos_kernels::impl::compute_seismograms( team_member.league_rank() * ChunkPolicy::tile_size * simd_size + tile; - if (starting_element_index >= nelements) { + if (starting_element_index >= nreceivers) { break; } const auto iterator = - policy.league_iterator(starting_element_index); + policy.mapped_league_iterator(starting_element_index); specfem::compute::load_on_device(team_member, iterator, field, element_field); diff --git a/include/policies/chunk.hpp b/include/policies/chunk.hpp index a4e53f93..9699d0df 100644 --- a/include/policies/chunk.hpp +++ b/include/policies/chunk.hpp @@ -52,6 +52,26 @@ struct chunk_index_type { const specfem::point::index index) : ielement(ielement), index(index){}; }; + +/** + * @brief Struct to store the index of a quadrature point generated by chunk + * policy. + * + * @tparam UseSIMD Indicates whether SIMD is used or not. + * @tparam DimensionType Dimension type of the elements within this iterator. + */ +template +struct mapped_chunk_index_type + : public chunk_index_type { + using Base = chunk_index_type; + int imap; ///< Index of the mapped element + + mapped_chunk_index_type(const int ielement, + const specfem::point::index index, + const int imap) + : Base(ielement, index), imap(imap) {} +}; + } // namespace impl /** @@ -93,7 +113,7 @@ class chunk { ///< type ///@} -private: +protected: constexpr static bool using_simd = simd::using_simd; constexpr static int simd_size = simd::size(); @@ -223,6 +243,35 @@ class chunk { return Kokkos::make_pair(indices(0), indices(num_elements - 1) + 1); } }; + +template +class mapped_chunk; + +template +class mapped_chunk + : public chunk { + using Base = chunk; + using mapped_index_type = + typename impl::mapped_chunk_index_type; ///< Index + +public: + mapped_chunk(const ViewType &indices, const ViewType &mapping, + const int ngllz, const int ngllx) + : Base(indices, ngllz, ngllx), mapping(mapping) {} + + KOKKOS_INLINE_FUNCTION + mapped_index_type operator()(const int i) const { + const auto base_index = Base::operator()(i); + return mapped_index_type(base_index.ielement, base_index.index, + mapping(base_index.ielement)); + } + +private: + ViewType mapping; +}; + } // namespace iterator namespace policy { @@ -237,7 +286,7 @@ template struct element_chunk : public Kokkos::TeamPolicy { -private: +protected: using IndexViewType = Kokkos::View< int *, typename ParallelConfig::execution_space::memory_space>; ///< View @@ -289,7 +338,7 @@ struct element_chunk true; ///< Indicates that this is a Kokkos team policy ///@} -private: +protected: constexpr static int simd_size = simd::size(); constexpr static bool using_simd = simd::using_simd; @@ -345,10 +394,44 @@ struct element_chunk return iterator_type(my_indices, ngllz, ngllx); } -private: +protected: IndexViewType elements; ///< View of elements int ngllz; ///< Number of GLL points in the z-direction int ngllx; ///< Number of GLL points in the x-direction }; + +template +struct mapped_element_chunk : public element_chunk { + using simd = typename ParallelConfig::simd; + using Base = element_chunk; + using IndexViewType = typename Base::IndexViewType; + + using mapped_iterator_type = + specfem::iterator::mapped_chunk; ///< Iterator + + mapped_element_chunk(const IndexViewType &view, const IndexViewType &mapping, + int ngllz, int ngllx) + : Base(view, ngllz, ngllx), mapping(mapping) {} + + KOKKOS_INLINE_FUNCTION + mapped_iterator_type mapped_league_iterator(const int start_index) const { + const int start = start_index; + const int end = + (start + Base::chunk_size * Base::simd_size > Base::elements.extent(0)) + ? Base::elements.extent(0) + : start + Base::chunk_size * Base::simd_size; + const auto elem_indices = + Kokkos::subview(Base::elements, Kokkos::make_pair(start, end)); + const auto map_indices = + Kokkos::subview(mapping, Kokkos::make_pair(start, end)); + return mapped_iterator_type(elem_indices, map_indices, Base::ngllz, + Base::ngllx); + } + +protected: + IndexViewType mapping; ///< View of elements +}; + } // namespace policy } // namespace specfem diff --git a/src/compute/compute_receivers.cpp b/src/compute/compute_receivers.cpp index 4e770a47..158af85e 100644 --- a/src/compute/compute_receivers.cpp +++ b/src/compute/compute_receivers.cpp @@ -24,19 +24,11 @@ specfem::compute::receivers::receivers( h_lagrange_interpolant(Kokkos::create_mirror_view(lagrange_interpolant)), elements("specfem::compute::receivers::elements", receivers.size()), h_elements(Kokkos::create_mirror_view(elements)), - receiver_domain_index_mapping( - "specfem::compute::receivers::receiver_domain_index_mapping", nspec), - h_receiver_domain_index_mapping( - Kokkos::create_mirror_view(receiver_domain_index_mapping)), element_types(element_types), impl::StationIterator(receivers.size(), stypes.size()), impl::SeismogramIterator(receivers.size(), stypes.size(), max_sig_step, dt, t0, nsteps_between_samples) { - for (int ispec = 0; ispec < nspec; ++ispec) { - h_receiver_domain_index_mapping(ispec) = -1; - } - for (int isies = 0; isies < stypes.size(); ++isies) { auto seis_type = stypes[isies]; switch (seis_type) { @@ -75,12 +67,6 @@ specfem::compute::receivers::receivers( }; const auto lcoord = specfem::algorithms::locate_point(gcoord, mesh); - if (h_receiver_domain_index_mapping(lcoord.ispec) != -1) { - throw std::runtime_error( - "Multiple receivers are detected in the same element"); - } - - h_receiver_domain_index_mapping(lcoord.ispec) = ireceiver; h_elements(ireceiver) = lcoord.ispec; const auto xi = mesh.quadratures.gll.h_xi; @@ -143,6 +129,17 @@ specfem::compute::receivers::receivers( GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)) = \ Kokkos::create_mirror_view( \ CREATE_VARIABLE_NAME(elements, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG))); \ + CREATE_VARIABLE_NAME(receiver_indices, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)) = \ + IndexViewType("specfem::compute::receivers::elements", \ + CREATE_VARIABLE_NAME(count, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), \ + GET_NAME(PROPERTY_TAG))); \ + CREATE_VARIABLE_NAME(h_receiver_indices, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)) = \ + Kokkos::create_mirror_view( \ + CREATE_VARIABLE_NAME(receiver_indices, GET_NAME(DIMENTION_TAG), \ GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG))); CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS( @@ -165,6 +162,11 @@ specfem::compute::receivers::receivers( (CREATE_VARIABLE_NAME(index, GET_NAME(DIMENTION_TAG), \ GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG))) = \ ispec; \ + CREATE_VARIABLE_NAME(h_receiver_indices, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)) \ + (CREATE_VARIABLE_NAME(index, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG))) = \ + ireceiver; \ CREATE_VARIABLE_NAME(index, GET_NAME(DIMENTION_TAG), \ GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)) \ ++; \ @@ -185,22 +187,24 @@ specfem::compute::receivers::receivers( Kokkos::deep_copy(lagrange_interpolant, h_lagrange_interpolant); Kokkos::deep_copy(elements, h_elements); - Kokkos::deep_copy(receiver_domain_index_mapping, - h_receiver_domain_index_mapping); return; } -Kokkos::View -specfem::compute::receivers::get_elements_on_host( +std::tuple, + Kokkos::View > +specfem::compute::receivers::get_indices_on_host( const specfem::element::medium_tag medium_tag, const specfem::element::property_tag property_tag) const { #define RETURN_VALUE(DIMENTION_TAG, MEDIUM_TAG, PROPERTY_TAG) \ if (medium_tag == GET_TAG(MEDIUM_TAG) && \ property_tag == GET_TAG(PROPERTY_TAG)) { \ - return CREATE_VARIABLE_NAME(h_elements, GET_NAME(DIMENTION_TAG), \ - GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)); \ + return std::make_tuple( \ + CREATE_VARIABLE_NAME(h_elements, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)), \ + CREATE_VARIABLE_NAME(h_receiver_indices, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG))); \ } CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS( @@ -211,16 +215,20 @@ specfem::compute::receivers::get_elements_on_host( #undef RETURN_VALUE } -Kokkos::View -specfem::compute::receivers::get_elements_on_device( +std::tuple, + Kokkos::View > +specfem::compute::receivers::get_indices_on_device( const specfem::element::medium_tag medium_tag, const specfem::element::property_tag property_tag) const { #define RETURN_VALUE(DIMENTION_TAG, MEDIUM_TAG, PROPERTY_TAG) \ if (medium_tag == GET_TAG(MEDIUM_TAG) && \ property_tag == GET_TAG(PROPERTY_TAG)) { \ - return CREATE_VARIABLE_NAME(elements, GET_NAME(DIMENTION_TAG), \ - GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)); \ + return std::make_tuple( \ + CREATE_VARIABLE_NAME(elements, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG)), \ + CREATE_VARIABLE_NAME(receiver_indices, GET_NAME(DIMENTION_TAG), \ + GET_NAME(MEDIUM_TAG), GET_NAME(PROPERTY_TAG))); \ } CALL_MACRO_FOR_ALL_MATERIAL_SYSTEMS( diff --git a/tests/unit-tests/displacement_tests/Newmark/serial/test1/STATIONS b/tests/unit-tests/displacement_tests/Newmark/serial/test1/STATIONS index 7f148ae7..f382abcb 100644 --- a/tests/unit-tests/displacement_tests/Newmark/serial/test1/STATIONS +++ b/tests/unit-tests/displacement_tests/Newmark/serial/test1/STATIONS @@ -3,6 +3,7 @@ S0002 AA 2250.0000000 3000.0000000 0.0 0.0 S0003 AA 2500.0000000 3000.0000000 0.0 0.0 S0004 AA 2750.0000000 3000.0000000 0.0 0.0 S0005 AA 3000.0000000 3000.0000000 0.0 0.0 +S0006 AA 2500.0000000 2500.0000000 0.0 0.0 S0007 AA 2500.0000000 2300.0000000 0.0 0.0 S0008 AA 2500.0000000 2100.0000000 0.0 0.0 S0009 AA 2500.0000000 1900.0000000 0.0 0.0