From 0692e7dc5a339dc205944587c8b6222e0ee84bce Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Sun, 18 Jan 2026 15:57:26 -0500 Subject: [PATCH 01/15] coupling integral wip -- pre-devel-merge --- .../dim2/impl/interface_container.hpp | 19 ++ .../dim2/nonconforming_interfaces.hpp | 37 +++- .../integrate/coupling_integral1d_dnshape.hpp | 166 ++++++++++++++++++ tests/unit-tests/CMakeLists.txt | 21 +++ .../dim2/coupling_integral/dshape.cpp | 105 +++++++++++ 5 files changed, 345 insertions(+), 3 deletions(-) create mode 100644 include/algorithms/integrate/coupling_integral1d_dnshape.hpp create mode 100644 tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp diff --git a/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp b/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp index df66a9087..dff149ebd 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp @@ -98,6 +98,25 @@ struct interface_container KOKKOS_FORCEINLINE_FUNCTION auto diff --git a/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp b/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp index 19ed53abe..c42b60e0b 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp @@ -66,9 +66,40 @@ class nonconforming_interfaces #ifndef NDEBUG // Debug check: abort if no matching specialization found - KOKKOS_ABORT_WITH_LOCATION( - "specfem::assembly::nonconforming_interfaces::get_interface_container(): No " - "matching specialization found."); + KOKKOS_ABORT_WITH_LOCATION("specfem::assembly::nonconforming_interfaces::" + "get_interface_container(): No " + "matching specialization found."); +#endif + + // Unreachable code - satisfy compiler return requirements + + SUPPRESS_TEMPORARY_REF(return {};) + } + + template + KOKKOS_INLINE_FUNCTION + InterfaceContainerType & + get_interface_container() { + // Compile-time dispatch using FOR_EACH_IN_PRODUCT macro + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), CONNECTION_TAG(NONCONFORMING), + INTERFACE_TAG(ELASTIC_ACOUSTIC, ACOUSTIC_ELASTIC), + BOUNDARY_TAG(NONE, ACOUSTIC_FREE_SURFACE, STACEY, + COMPOSITE_STACEY_DIRICHLET)), + CAPTURE((interface_container, interface_container)) { + if constexpr (InterfaceTag == _interface_tag_ && + BoundaryTag == _boundary_tag_ && + ConnectionTag == _connection_tag_) { + return _interface_container_; + } + }) + +#ifndef NDEBUG + // Debug check: abort if no matching specialization found + KOKKOS_ABORT_WITH_LOCATION("specfem::assembly::nonconforming_interfaces::" + "get_interface_container(): No " + "matching specialization found."); #endif // Unreachable code - satisfy compiler return requirements diff --git a/include/algorithms/integrate/coupling_integral1d_dnshape.hpp b/include/algorithms/integrate/coupling_integral1d_dnshape.hpp new file mode 100644 index 000000000..3fe90bd8e --- /dev/null +++ b/include/algorithms/integrate/coupling_integral1d_dnshape.hpp @@ -0,0 +1,166 @@ +#pragma once + +#include "enumerations/interface.hpp" +#include "execution/for_each_level.hpp" +#include "execution/team_thread_md_range_iterator.hpp" +#include "specfem/assembly.hpp" +#include "specfem/point.hpp" +#include + +/** + * @brief Temporary file for integral in issue # 1550. this function should be + * merged with `coupling_integral1d.hpp:coupling_integral` once the intersection + * field typing is updated to store two fields. + * + */ + +namespace specfem::algorithms { + +/** + * @brief Takes a field `intersection_field` on the intersection and computes, + * for each self GLL point, the integral of `intersection_field` times the + * normal derivative of the shape function at that point. `intersection_field` + * should be call-accessible (e.g. Kokkos::View) with shape: + * + * (chunk_size, n_quad_intersection, self::components()) + * + * After handling any other intersection forces, boundary conditions, etc. the + * result can be `atomic_add`ed to the acceleration field. + * + * @tparam dimension_tag dimension of the simulation + * @tparam IndexType The chunk_edge iterator type + * @tparam IntersectionFieldViewType The type of `intersection_field` + * @tparam ChunkEdgeWeightJacobianType A nonconforming chunk_edge accessor + * holding `intersection_factor` + * @tparam CallableType The callback function, which will be given the point + * index and corresponding evaluated integral + * @param nonconforming_interfaces - assembly.nonconforming_interfaces struct + * @param lagrange_derivative - Local derivatives of shape functions on the + * element + * @param chunk_index - the outer index (chunk_edge) that gets iterated for + * points + * @param intersection_field - the field to integrate + * @param weight_jacobian - nonconforming chunk_edge accessor holding + * `intersection_factor` + * @param intersection_normal_covariant_edgelocal - the normal vector at each + * intersection point in contravariant local coordinate basis, with the first + * component in the normal direction. + * @param callback - callback function to capture integral values + * @ingroup AlgorithmsIntegration + */ +template +KOKKOS_FUNCTION void coupling_integral_dnshape( + const specfem::assembly::nonconforming_interfaces + &nonconforming_interfaces, + const QuadratureType &lagrange_derivative, const IndexType &chunk_index, + const IntersectionFieldViewType &intersection_field, + const IntersectionFactor &intersection_factor, + const IntersectionJacobianType &intersection_normal_covariant_edgelocal, + const TransferFunctionDerivativeType + &transfer_function_self_derivative /* This could be loaded just like + transfer_function_self */ + , + const CallableType &callback) { + + constexpr auto self_medium_tag = specfem::interface::attributes< + dimension_tag, IntersectionFactor::interface_tag>::self_medium(); + + using PointIndexType = + typename IndexType::iterator_type::index_type::index_type; + using PointFieldType = + specfem::point::acceleration; + using SelfTransferFunctionType = specfem::point::transfer_function_self< + IntersectionFactor::n_quad_intersection, dimension_tag, + IntersectionFactor::interface_tag, IntersectionFactor::boundary_tag>; + + constexpr int ncomp = PointFieldType::components; + constexpr int nquad_intersection = IntersectionFactor::n_quad_intersection; + constexpr int nquad_element = IntersectionFactor::n_quad_element; + + specfem::execution::for_each_level( + chunk_index.get_iterator(), + [&](const typename IndexType::iterator_type::index_type &index) { + const auto self_index = index.get_index(); + const auto self_index_local = index.get_local_index(); + const bool ipoint_n_is_ix = + self_index.edge_type == specfem::mesh_entity::dim2::type::left || + self_index.edge_type == specfem::mesh_entity::dim2::type::right; + const int &iedge = self_index_local.iedge; + SelfTransferFunctionType transfer_function_self; + specfem::assembly::load_on_device(self_index, nonconforming_interfaces, + transfer_function_self); + + const int &ipoint_s = self_index.ipoint; + for (int ipoint_n = 0; ipoint_n < nquad_element; ++ipoint_n) { + // iterate backwards from the edge (this gets called for each gllxz) + specfem::point::index interior_index( + self_index.ispec, self_index.iz, self_index.ix); + + // sample derivative of interior_index shape function, in normal + // (covector) direction, at edge coordinate. + const type_real dlagrange_dn = + ipoint_n_is_ix + ? lagrange_derivative.xi(ipoint_n, self_index.ix) + : lagrange_derivative.gamma(ipoint_n, self_index.iz); + + // we may want to refactor this at some point. For now, this should be + // fine to update the index point. + switch (self_index.edge_type) { + case specfem::mesh_entity::dim2::type::right: + interior_index.ix = (nquad_element - 1) - ipoint_n; + break; + case specfem::mesh_entity::dim2::type::top: + interior_index.iz = (nquad_element - 1) - ipoint_n; + break; + case specfem::mesh_entity::dim2::type::left: + interior_index.ix = ipoint_n; + break; + case specfem::mesh_entity::dim2::type::bottom: + interior_index.iz = ipoint_n; + break; + } + + PointFieldType result; +#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) +#pragma unroll +#endif + for (int icomp = 0; icomp < ncomp; icomp++) { + result(icomp) = 0; + } +#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) +#pragma unroll +#endif + for (int iquad = 0; iquad < nquad_intersection; iquad++) { + +#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) +#pragma unroll +#endif + for (int icomp = 0; icomp < ncomp; icomp++) { + // dshape_dn, local-normal derivative (derivative of + // normal-direction L, which is normally kronecker delta + // indicating on edge) + type_real dshape_dn = + intersection_normal_covariant_edgelocal(iedge, iquad, 0) * + (dlagrange_dn * transfer_function_self(iquad)); + // dshape_dn, local-tangential derivative (differentiate + // transfer_function_self instead of normal-direction L) + if (ipoint_n == 0) { + dshape_dn += + intersection_normal_covariant_edgelocal(iedge, iquad, 1) * + transfer_function_self_derivative(iedge, ipoint_s, iquad); + } + result(icomp) += intersection_field(iedge, iquad, icomp) * + intersection_factor(iedge, iquad) * dshape_dn; + } + } + + callback(self_index, result); + } + }); +} + +} // namespace specfem::algorithms diff --git a/tests/unit-tests/CMakeLists.txt b/tests/unit-tests/CMakeLists.txt index 13bdb46bb..f04cacba0 100644 --- a/tests/unit-tests/CMakeLists.txt +++ b/tests/unit-tests/CMakeLists.txt @@ -626,6 +626,27 @@ target_link_libraries( -lpthread -lm ) +add_executable( + coupling_integral_tests + algorithms/dim2/coupling_integral/dshape.cpp +) + +target_link_libraries( + coupling_integral_tests + algorithms + enumerations + mesh + assembly + gtest_main + Kokkos::kokkos + specfem_environment + utilities + + quadrature + yaml-cpp + -lpthread -lm +) + add_executable( policies_tests policies/policies.cpp diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp new file mode 100644 index 000000000..881c42242 --- /dev/null +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -0,0 +1,105 @@ +#include "algorithms/integrate/coupling_integral1d_dnshape.hpp" + +#include "SPECFEM_Environment.hpp" +#include + +class ChunkEdgeIndex { +public: + static constexpr auto accessor_type = + specfem::data_access::AccessorType::chunk_edge; + using KokkosIndexType = Kokkos::TeamPolicy<>::member_type; + + /** + * @brief Get Kokkos team member index. + * @return Reference to Kokkos team member + */ + KOKKOS_INLINE_FUNCTION + constexpr const KokkosIndexType &get_policy_index() const { + return this->kokkos_index; + } + + /** + * @brief Construct chunk edge index. + * @param nedges Number of edges in chunk + * @param kokkos_index Kokkos team member + */ + KOKKOS_INLINE_FUNCTION + ChunkEdgeIndex(const int nedges, const KokkosIndexType &kokkos_index) + : kokkos_index(kokkos_index), _nedges(nedges) {} + + /** + * @brief Get number of edges. + * @return Edge count + */ + KOKKOS_INLINE_FUNCTION int nedges() const { return _nedges; } + +private: + int _nedges; ///< Number of edges in the chunk + KokkosIndexType kokkos_index; /**< Kokkos team member for this chunk */ +}; + +// temporary test for purposes of uncombined coupling_integral +TEST(CouplingIntegral, SimpleDShapeTest) { + constexpr int nquad_element = 5; + constexpr int nquad_intersection = 5; + constexpr int ngllx = nquad_element; + constexpr int ngllz = nquad_element; + constexpr int num_edges = 1; + + constexpr auto interface_tag = + specfem::interface::interface_tag::acoustic_elastic; + constexpr auto boundary_tag = specfem::element::boundary_tag::none; + using memory_space = Kokkos::DefaultExecutionSpace::memory_space; + + constexpr auto medium_self = + specfem::interface::attributes::self_medium(); + constexpr auto ncomp_self = + specfem::element::attributes::components; + + specfem::assembly::nonconforming_interfaces + nonconforming_interfaces; + + nonconforming_interfaces.template get_interface_container< + interface_tag, boundary_tag, + specfem::connections::type::nonconforming>() = + specfem::assembly::nonconforming_interfaces_impl::interface_container< + specfem::dimension::type::dim2, interface_tag, boundary_tag, + specfem::connections::type::nonconforming>(ngllx, ngllz, num_edges); + + using TransferFunctionType = specfem::chunk_edge::impl::transfer_function< + dimension_tag, 1, nquad_intersection, nquad_edge, + specfem::data_access::DataClassType::transfer_function_self, + interface_tag, boundary_tag, memory_space, Kokkos::MemoryTraits<> >; + using FunctionType = specfem::datatype::VectorChunkEdgeViewType< + type_real, dimension_tag, 1, nquad_intersection, ncomp_self, false, + memory_space, Kokkos::MemoryTraits<> >; + + Kokkos::parallel_for( + "SimpleDShapeTest", Kokkos::TeamPolicy<>(num_edges, 1, 1), + KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type &team_member) { + const int iedge = team_member.league_rank(); + const TransferFunctionType TF(Kokkos::subview( + transfer_function_view, Kokkos::make_pair(iedge, iedge + 1), + Kokkos::ALL(), Kokkos::ALL())); + const FunctionType F( + Kokkos::subview(function_view, Kokkos::make_pair(iedge, iedge + 1), + Kokkos::ALL(), Kokkos::ALL())); + specfem::algorithms::coupling_integral_dnshape( + nonconforming_interfaces, ChunkEdgeIndex(num_edges, team_member), + TF, F, [&](const auto &index, const auto &point) { + for (int icomp = 0; icomp < ncomp_self; ++icomp) { + Kokkos::single(Kokkos::PerTeam(team_member), [&]() { + result_view(index(0), index(1), icomp) = point(icomp); + }); + } + }); + }); +} + +int main(int argc, char *argv[]) { + ::testing::InitGoogleTest(&argc, argv); + ::testing::AddGlobalTestEnvironment(new SPECFEMEnvironment); + return RUN_ALL_TESTS(); +} From cd215a017e6fa4a8e935bd3fc7390bdcca9f2430 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Sun, 18 Jan 2026 23:35:17 -0500 Subject: [PATCH 02/15] kernel up to validation --- .../dim2/impl/interface_container.hpp | 11 +- .../dim2/nonconforming_interfaces.hpp | 33 +- .../integrate/coupling_integral1d_dnshape.hpp | 27 +- .../dim2/coupling_integral/dshape.cpp | 332 +++++++++++++++--- .../include/fixture/impl/accessors.hpp | 44 ++- .../fixture/nonconforming_interface.hpp | 1 + .../analytical_function.hpp | 102 ++++++ .../nonconforming_interface/initializers.hpp | 10 + .../lagrange_derivative.hpp | 49 +++ .../nonconforming_interface/quadrature.hpp | 148 ++++++-- 10 files changed, 639 insertions(+), 118 deletions(-) create mode 100644 tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp diff --git a/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp b/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp index dff149ebd..c5036eb8e 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp @@ -98,22 +98,23 @@ struct interface_container public: static constexpr auto dimension_tag = specfem::dimension::type::dim2; -private: +protected: template @@ -64,37 +64,6 @@ class nonconforming_interfaces } }) -#ifndef NDEBUG - // Debug check: abort if no matching specialization found - KOKKOS_ABORT_WITH_LOCATION("specfem::assembly::nonconforming_interfaces::" - "get_interface_container(): No " - "matching specialization found."); -#endif - - // Unreachable code - satisfy compiler return requirements - - SUPPRESS_TEMPORARY_REF(return {};) - } - - template - KOKKOS_INLINE_FUNCTION - InterfaceContainerType & - get_interface_container() { - // Compile-time dispatch using FOR_EACH_IN_PRODUCT macro - FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), CONNECTION_TAG(NONCONFORMING), - INTERFACE_TAG(ELASTIC_ACOUSTIC, ACOUSTIC_ELASTIC), - BOUNDARY_TAG(NONE, ACOUSTIC_FREE_SURFACE, STACEY, - COMPOSITE_STACEY_DIRICHLET)), - CAPTURE((interface_container, interface_container)) { - if constexpr (InterfaceTag == _interface_tag_ && - BoundaryTag == _boundary_tag_ && - ConnectionTag == _connection_tag_) { - return _interface_container_; - } - }) - #ifndef NDEBUG // Debug check: abort if no matching specialization found KOKKOS_ABORT_WITH_LOCATION("specfem::assembly::nonconforming_interfaces::" diff --git a/include/algorithms/integrate/coupling_integral1d_dnshape.hpp b/include/algorithms/integrate/coupling_integral1d_dnshape.hpp index 3fe90bd8e..00cd1f1de 100644 --- a/include/algorithms/integrate/coupling_integral1d_dnshape.hpp +++ b/include/algorithms/integrate/coupling_integral1d_dnshape.hpp @@ -42,9 +42,14 @@ namespace specfem::algorithms { * @param intersection_field - the field to integrate * @param weight_jacobian - nonconforming chunk_edge accessor holding * `intersection_factor` - * @param intersection_normal_covariant_edgelocal - the normal vector at each - * intersection point in contravariant local coordinate basis, with the first - * component in the normal direction. + * @param intersection_normal_contravariant_edgelocal - the normal vector at + * each intersection point in contravariant local coordinate basis, with the + * first component in the outward normal direction (-/+xi for left/right, + * +/-gamma for top/bottom), and second component in the tangential axis (+gamma + * for left/right, +xi for top/bottom). + * @param transfer_function_self_derivative - (TEMPORARY, TO BE LOADED LATER BY + * THIS FUNCTION) derivative (in edge coordinates, not intersection coordinates + * ) of transfer_function_self, evaluated at intersection quadrature points. * @param callback - callback function to capture integral values * @ingroup AlgorithmsIntegration */ @@ -55,10 +60,11 @@ template &nonconforming_interfaces, + const int &ngllz, const int &ngllx, const QuadratureType &lagrange_derivative, const IndexType &chunk_index, const IntersectionFieldViewType &intersection_field, const IntersectionFactor &intersection_factor, - const IntersectionJacobianType &intersection_normal_covariant_edgelocal, + const IntersectionJacobianType &intersection_normal_contravariant_edgelocal, const TransferFunctionDerivativeType &transfer_function_self_derivative /* This could be loaded just like transfer_function_self */ @@ -79,7 +85,6 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( constexpr int ncomp = PointFieldType::components; constexpr int nquad_intersection = IntersectionFactor::n_quad_intersection; - constexpr int nquad_element = IntersectionFactor::n_quad_element; specfem::execution::for_each_level( chunk_index.get_iterator(), @@ -95,7 +100,8 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( transfer_function_self); const int &ipoint_s = self_index.ipoint; - for (int ipoint_n = 0; ipoint_n < nquad_element; ++ipoint_n) { + const int ngll_n = ipoint_n_is_ix ? ngllx : ngllz; + for (int ipoint_n = 0; ipoint_n < ngll_n; ++ipoint_n) { // iterate backwards from the edge (this gets called for each gllxz) specfem::point::index interior_index( self_index.ispec, self_index.iz, self_index.ix); @@ -111,10 +117,10 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( // fine to update the index point. switch (self_index.edge_type) { case specfem::mesh_entity::dim2::type::right: - interior_index.ix = (nquad_element - 1) - ipoint_n; + interior_index.ix = (ngllx - 1) - ipoint_n; break; case specfem::mesh_entity::dim2::type::top: - interior_index.iz = (nquad_element - 1) - ipoint_n; + interior_index.iz = (ngllz - 1) - ipoint_n; break; case specfem::mesh_entity::dim2::type::left: interior_index.ix = ipoint_n; @@ -144,13 +150,14 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( // normal-direction L, which is normally kronecker delta // indicating on edge) type_real dshape_dn = - intersection_normal_covariant_edgelocal(iedge, iquad, 0) * + intersection_normal_contravariant_edgelocal(iedge, iquad, 0) * (dlagrange_dn * transfer_function_self(iquad)); // dshape_dn, local-tangential derivative (differentiate // transfer_function_self instead of normal-direction L) if (ipoint_n == 0) { dshape_dn += - intersection_normal_covariant_edgelocal(iedge, iquad, 1) * + intersection_normal_contravariant_edgelocal(iedge, iquad, + 1) * transfer_function_self_derivative(iedge, ipoint_s, iquad); } result(icomp) += intersection_field(iedge, iquad, icomp) * diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index 881c42242..a8ba63b09 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -1,13 +1,83 @@ #include "algorithms/integrate/coupling_integral1d_dnshape.hpp" +#include "specfem/chunk_edge.hpp" + +#include "utilities/include/fixture/impl/accessors.hpp" +#include "utilities/include/fixture/nonconforming_interface.hpp" #include "SPECFEM_Environment.hpp" +#include "utilities/include/fixture/nonconforming_interface/analytical_function.hpp" +#include "utilities/include/fixture/nonconforming_interface/intersection_function.hpp" +#include "utilities/include/fixture/nonconforming_interface/quadrature.hpp" +#include "utilities/include/fixture/nonconforming_interface/transfer_function.hpp" #include -class ChunkEdgeIndex { +template +class ChunkEdgeIterator : public specfem::execution::TeamThreadRangePolicy< + Kokkos::TeamPolicy<>::member_type, int> { +private: + using KokkosIndexType = Kokkos::TeamPolicy<>::member_type; + +public: + using base_type = + specfem::execution::TeamThreadRangePolicy; + using execution_space = typename base_type::execution_space; + using index_type = + specfem::execution::EdgePointIndex; + + KOKKOS_INLINE_FUNCTION + ChunkEdgeIterator(const EdgeTypesView &edge_types, const int &nedges, + const int &ngllz, const int &ngllx, + const KokkosIndexType &team_member) + : edge_types(edge_types), nedges(nedges), ngllz(ngllz), ngllx(ngllz), + num_points(std::max(ngllx, ngllz)), + base_type(team_member, nedges * std::max(ngllx, ngllz)), + global_offset(0) {} + + KOKKOS_INLINE_FUNCTION + const index_type + operator()(const typename base_type::policy_index_type &i) const { + const int iedge = i % nedges; + const int ipoint = i / nedges; + + const specfem::mesh_entity::dim2::type edge_type = edge_types(iedge); + const bool is_leftright = + (edge_type == specfem::mesh_entity::dim2::type::left || + edge_type == specfem::mesh_entity::dim2::type::right); + const int num_points_norm = is_leftright ? ngllx : ngllz; + + const int inorm = (edge_type == specfem::mesh_entity::dim2::type::bottom || + edge_type == specfem::mesh_entity::dim2::type::left) + ? 0 + : num_points_norm - 1; + + const int iz = is_leftright ? ipoint : inorm; + const int ix = is_leftright ? inorm : ipoint; + + return index_type( + specfem::point::edge_index( + 0 /*ispec*/, global_offset + iedge /*iedge*/, ipoint, iz, ix, + edge_type), + iedge, i); + } + + const int nedges; + +private: + EdgeTypesView edge_types; + const int global_offset; + const int ngllz; + const int ngllx; + const int num_points; +}; + +template class ChunkEdgeIndex { public: static constexpr auto accessor_type = specfem::data_access::AccessorType::chunk_edge; using KokkosIndexType = Kokkos::TeamPolicy<>::member_type; + using iterator_type = ChunkEdgeIterator; /** * @brief Get Kokkos team member index. @@ -24,8 +94,11 @@ class ChunkEdgeIndex { * @param kokkos_index Kokkos team member */ KOKKOS_INLINE_FUNCTION - ChunkEdgeIndex(const int nedges, const KokkosIndexType &kokkos_index) - : kokkos_index(kokkos_index), _nedges(nedges) {} + ChunkEdgeIndex(const EdgeTypesView &edge_types, const int &nedges, + const int &ngllz, const int &ngllx, + const KokkosIndexType &kokkos_index) + : kokkos_index(kokkos_index), _nedges(nedges), + iterator(edge_types, nedges, ngllz, ngllx, kokkos_index) {} /** * @brief Get number of edges. @@ -36,64 +109,241 @@ class ChunkEdgeIndex { private: int _nedges; ///< Number of edges in the chunk KokkosIndexType kokkos_index; /**< Kokkos team member for this chunk */ + iterator_type iterator; + +public: + KOKKOS_INLINE_FUNCTION const iterator_type &get_iterator() const { + return iterator; + } +}; + +/** + * @brief Patches assembly::nonconforming_interfaces to not require mesh and + * edge types. + * + * nonconforming_interfaces only provides const access to the impl containers. + * This class grants access to resizing these containers directly. Note that + * const access still allows modification of values inside the views. + */ +class nonconforming_interfaces_patch + : public specfem::assembly::nonconforming_interfaces< + specfem::dimension::type::dim2> { + int ngllz; + int ngllx; + int nquad_intersection; + +public: + nonconforming_interfaces_patch(const int &ngllz, const int &ngllx, + const int &nquad_intersection) + : ngllz(ngllz), ngllx(ngllx), nquad_intersection(nquad_intersection) {}; + + template + void reinit_container(const int &num_edges) { + + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM2), CONNECTION_TAG(NONCONFORMING), + INTERFACE_TAG(ELASTIC_ACOUSTIC, ACOUSTIC_ELASTIC), + BOUNDARY_TAG(NONE, STACEY, ACOUSTIC_FREE_SURFACE, + COMPOSITE_STACEY_DIRICHLET)), + CAPTURE(interface_container) { + if constexpr (_interface_tag_ == InterfaceTag && + _boundary_tag_ == BoundaryTag && + _connection_tag_ == ConnectionTag) { + _interface_container_ = + InterfaceContainerType<_interface_tag_, _boundary_tag_, + _connection_tag_>( + ngllz, ngllx, nquad_intersection, num_edges); + } + }) + } }; // temporary test for purposes of uncombined coupling_integral TEST(CouplingIntegral, SimpleDShapeTest) { - constexpr int nquad_element = 5; - constexpr int nquad_intersection = 5; - constexpr int ngllx = nquad_element; - constexpr int ngllz = nquad_element; - constexpr int num_edges = 1; - + constexpr auto dimension_tag = specfem::dimension::type::dim2; constexpr auto interface_tag = specfem::interface::interface_tag::acoustic_elastic; constexpr auto boundary_tag = specfem::element::boundary_tag::none; using memory_space = Kokkos::DefaultExecutionSpace::memory_space; + // ========================================================== + // test case (this will be replaced by fixtures eventually) + + using QuadX = specfem::test_fixture::QuadraturePoints::GLL2; + using QuadZ = specfem::test_fixture::QuadraturePoints::GLL2; + using QuadIntersection = specfem::test_fixture::QuadraturePoints::Asymm5Point; + constexpr int num_edges = 1; + + const type_real intersection_min = -1; + const type_real intersection_max = 1; + const auto side = specfem::mesh_entity::dim2::type::top; + // we are integrating d(shape-function)/dn * F, where + // :: F(x) = x^2 + // :: n = [1, x] in the contravariant local coordinate basis (see + // coupling_integral1d_dnshape.hpp) + // :: shape functions on GLL2 x GLL1 elements + // :: integral on [-1,1] (intersection_) + // :: edge has same coordinates, and goes between [-1,1] in xi + + // structs to use + + using LagrangeDerivative2D = + specfem::test_fixture::LagrangeDerivative2D; + LagrangeDerivative2D lagrange_derivative("dshape::lagrange_derivative"); + + using IntersectionFunctionInitializer = specfem::test_fixture:: + IntersectionFunctionInitializer2D::FromAnalyticalFunction< + specfem::test_fixture::AnalyticalFunctionType::Power<2>, + QuadIntersection>; + const auto function_view = specfem::test_fixture::IntersectionFunction2D( + IntersectionFunctionInitializer()) + .get_view(); + + // we will probably fixturize this at some point to make it not so bloated + using IntersectionFactor = + specfem::test_fixture::impl::NonconformingIntersectionFactorPatch< + interface_tag, boundary_tag, num_edges, QuadIntersection::nquad>; + auto intersection_factor = IntersectionFactor("dshape::intersection_factor"); + + { + const auto h_intersection_factor = intersection_factor.create_host_mirror(); + + for (int iedge = 0; iedge < num_edges; iedge++) { + for (int iquad = 0; iquad < QuadIntersection::nquad; iquad++) { + h_intersection_factor(iedge, iquad) = + specfem::test_fixture::QuadratureRule< + QuadIntersection>::compute_lagrange_quadrature_weight(iquad, -1, + 1); + } + } + intersection_factor.sync_to_device(h_intersection_factor); + } + + // Contravariant normal derivative is not a datatype yet. Keep this for now. + using IntersectionContraNormalFunction = + specfem::test_fixture::AnalyticalFunctionType::Chain< + specfem::test_fixture::AnalyticalFunctionType::Power<0>, + specfem::test_fixture::AnalyticalFunctionType::Power<1> >; + using IntersectionContraNormalInitializer = + specfem::test_fixture::IntersectionFunctionInitializer2D:: + FromAnalyticalFunction; + const auto intersection_contra_normal = + specfem::test_fixture::IntersectionFunction2D( + IntersectionContraNormalInitializer()) + .get_view(); + + // no data access help for transfer_function_self derivative, but the type + // only needs to support deriv(iedge, iquad, icomp). + using TransferFunctionDerivativeType = + specfem::datatype::VectorChunkEdgeViewType< + type_real, dimension_tag, num_edges, QuadIntersection::nquad, 2, + false, memory_space, Kokkos::MemoryTraits<> >; + + TransferFunctionDerivativeType transfer_function_self_derivative( + "dshape::transfer_function_self_derivative"); // init later with + // interface_container + + // ========================================================== + + constexpr int nquad_intersection = QuadIntersection::nquad; + constexpr int ngllx = QuadX::nquad; + constexpr int ngllz = QuadZ::nquad; + constexpr int nquad_element = std::max(ngllx, ngllz); + constexpr int chunk_size = 1; + constexpr auto medium_self = - specfem::interface::attributes::self_medium(); constexpr auto ncomp_self = - specfem::element::attributes::components; - - specfem::assembly::nonconforming_interfaces - nonconforming_interfaces; - - nonconforming_interfaces.template get_interface_container< - interface_tag, boundary_tag, - specfem::connections::type::nonconforming>() = - specfem::assembly::nonconforming_interfaces_impl::interface_container< - specfem::dimension::type::dim2, interface_tag, boundary_tag, - specfem::connections::type::nonconforming>(ngllx, ngllz, num_edges); - - using TransferFunctionType = specfem::chunk_edge::impl::transfer_function< - dimension_tag, 1, nquad_intersection, nquad_edge, - specfem::data_access::DataClassType::transfer_function_self, - interface_tag, boundary_tag, memory_space, Kokkos::MemoryTraits<> >; + specfem::element::attributes::components; + + nonconforming_interfaces_patch nonconforming_interfaces(ngllz, ngllx, + nquad_intersection); + nonconforming_interfaces.template reinit_container< + interface_tag, boundary_tag, specfem::connections::type::nonconforming>( + num_edges); + + const auto &interface_container = + nonconforming_interfaces.template get_interface_container< + interface_tag, boundary_tag, + specfem::connections::type::nonconforming>(); + + // ================================================================= + // populate this nonconforming interface container and + // transfer_function_self_derivative + { + TransferFunctionDerivativeType::HostMirror h_tfsd = + Kokkos::create_mirror_view(transfer_function_self_derivative); + for (int iedge = 0; iedge < num_edges; iedge++) { + for (int iquad_edge = 0; iquad_edge < ngllx; iquad_edge++) { + for (int iquad_intersection = 0; + iquad_intersection < nquad_intersection; iquad_intersection++) { + + // no transformation (it would go here, otherwise) + const double intersection_point_in_edge_coords = + QuadIntersection::quadrature_points[iquad_intersection]; + interface_container.h_transfer_function(iedge, iquad_intersection, + iquad_edge) = specfem:: + test_fixture::QuadratureRule::evaluate_lagrange_polynomial( + iquad_edge, intersection_point_in_edge_coords); + h_tfsd(iedge, iquad_edge, iquad_intersection) = specfem:: + test_fixture::QuadratureRule::evaluate_lagrange_derivative( + iquad_edge, intersection_point_in_edge_coords); + } + } + } + Kokkos::deep_copy(transfer_function_self_derivative, h_tfsd); + Kokkos::deep_copy(interface_container.transfer_function, + interface_container.h_transfer_function); + } + // ================================================================= + + using EdgeTypesView = Kokkos::View >; + EdgeTypesView edge_types("dshape::edge_types", num_edges); + EdgeTypesView::HostMirror h_edge_types = + Kokkos::create_mirror_view(edge_types); + for (int iedge = 0; iedge < num_edges; iedge++) { + h_edge_types(iedge) = side; + } + Kokkos::deep_copy(edge_types, h_edge_types); + + // chunk subviews (FunctionType comes from copied test template -- skip the + // rest and assume num_edges == 1) using FunctionType = specfem::datatype::VectorChunkEdgeViewType< - type_real, dimension_tag, 1, nquad_intersection, ncomp_self, false, - memory_space, Kokkos::MemoryTraits<> >; + type_real, dimension_tag, chunk_size, nquad_intersection, ncomp_self, + false, memory_space, Kokkos::MemoryTraits<> >; + + const int num_chunks = + num_edges / chunk_size + ((num_edges % chunk_size != 0) ? 1 : 0); Kokkos::parallel_for( "SimpleDShapeTest", Kokkos::TeamPolicy<>(num_edges, 1, 1), KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type &team_member) { const int iedge = team_member.league_rank(); - const TransferFunctionType TF(Kokkos::subview( - transfer_function_view, Kokkos::make_pair(iedge, iedge + 1), + const FunctionType F(Kokkos::subview( + function_view, + Kokkos::make_pair(iedge * chunk_size, (iedge + 1) * chunk_size), Kokkos::ALL(), Kokkos::ALL())); - const FunctionType F( - Kokkos::subview(function_view, Kokkos::make_pair(iedge, iedge + 1), - Kokkos::ALL(), Kokkos::ALL())); specfem::algorithms::coupling_integral_dnshape( - nonconforming_interfaces, ChunkEdgeIndex(num_edges, team_member), - TF, F, [&](const auto &index, const auto &point) { - for (int icomp = 0; icomp < ncomp_self; ++icomp) { - Kokkos::single(Kokkos::PerTeam(team_member), [&]() { - result_view(index(0), index(1), icomp) = point(icomp); - }); - } + nonconforming_interfaces, ngllz, ngllx, lagrange_derivative, + ChunkEdgeIndex( + Kokkos::subview(edge_types, + Kokkos::make_pair(iedge * chunk_size, + (iedge + 1) * chunk_size)), + num_edges, ngllz, ngllx, team_member), + F, intersection_factor, intersection_contra_normal, + transfer_function_self_derivative, + [&](const auto &index, const auto &point) { + // TODO START HERE + // for (int icomp = 0; icomp < ncomp_self; ++icomp) { + // Kokkos::single(Kokkos::PerTeam(team_member), [&]() { + // result_view(index(0), index(1), icomp) = point(icomp); + // }); + // } }); }); } diff --git a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp index 63b5af130..3959afce7 100644 --- a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp +++ b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp @@ -1,9 +1,11 @@ #pragma once +#include "Kokkos_Core.hpp" #include "enumerations/coupled_interface.hpp" #include "enumerations/medium.hpp" #include "specfem/data_access/accessor.hpp" #include "specfem/data_access/data_class.hpp" +#include namespace specfem::test_fixture::impl { /** @@ -29,8 +31,22 @@ struct NonconformingAccessorPatch2D static constexpr auto connection_tag = specfem::connections::type::nonconforming; /// View type for storing intersection scaling factors +private: + template struct ViewInternalType; + + template + struct ViewInternalType { + using type = + typename ViewInternalType::type; + }; + + template struct ViewInternalType { + using type = LeftType; + }; + +public: using DataViewType = - Kokkos::View::type, Kokkos::DefaultExecutionSpace::memory_space>; private: @@ -53,6 +69,14 @@ struct NonconformingAccessorPatch2D KOKKOS_INLINE_FUNCTION auto &operator()(Indices... indices) const { return data_(indices...); } + + typename DataViewType::HostMirror create_host_mirror() { + return Kokkos::create_mirror_view(data_); + } + + void sync_to_device(const typename DataViewType::HostMirror &host_data) { + Kokkos::deep_copy(data_, host_data); + } }; template { + using NonconformingAccessorPatch2D< + InterfaceTag, BoundaryTag, + specfem::data_access::DataClassType::transfer_function_self, + NumberElements, NQuadElement, + NQuadIntersection>::NonconformingAccessorPatch2D; static constexpr int chunk_size = NumberElements; static constexpr int n_quad_element = NQuadElement; static constexpr int n_quad_intersection = NQuadIntersection; @@ -75,6 +104,11 @@ struct NonconformingTransferFunctionCoupledPatch InterfaceTag, BoundaryTag, specfem::data_access::DataClassType::transfer_function_coupled, NumberElements, NQuadElement, NQuadIntersection> { + using NonconformingAccessorPatch2D< + InterfaceTag, BoundaryTag, + specfem::data_access::DataClassType::transfer_function_coupled, + NumberElements, NQuadElement, + NQuadIntersection>::NonconformingAccessorPatch2D; static constexpr int chunk_size = NumberElements; static constexpr int n_quad_element = NQuadElement; static constexpr int n_quad_intersection = NQuadIntersection; @@ -87,6 +121,10 @@ struct NonconformingIntersectionNormalPatch InterfaceTag, BoundaryTag, specfem::data_access::DataClassType::intersection_normal, NumberElements, NQuadIntersection, 2> { + using NonconformingAccessorPatch2D< + InterfaceTag, BoundaryTag, + specfem::data_access::DataClassType::intersection_normal, NumberElements, + NQuadIntersection, 2>::NonconformingAccessorPatch2D; static constexpr int chunk_size = NumberElements; static constexpr int n_quad_intersection = NQuadIntersection; }; @@ -98,6 +136,10 @@ struct NonconformingIntersectionFactorPatch InterfaceTag, BoundaryTag, specfem::data_access::DataClassType::intersection_factor, NumberElements, NQuadIntersection> { + using NonconformingAccessorPatch2D< + InterfaceTag, BoundaryTag, + specfem::data_access::DataClassType::intersection_factor, NumberElements, + NQuadIntersection>::NonconformingAccessorPatch2D; static constexpr int chunk_size = NumberElements; static constexpr int n_quad_intersection = NQuadIntersection; }; diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface.hpp index d59d61b8b..b836d7dc9 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface.hpp @@ -5,5 +5,6 @@ #include "nonconforming_interface/edge_function.hpp" #include "nonconforming_interface/intersection_data.hpp" #include "nonconforming_interface/intersection_function.hpp" +#include "nonconforming_interface/lagrange_derivative.hpp" #include "nonconforming_interface/quadrature.hpp" #include "nonconforming_interface/transfer_function.hpp" diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp index 699d7ce1e..837ed2a22 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp @@ -1,5 +1,7 @@ #pragma once +#include "../impl/descriptions.hpp" +#include "Serial/Kokkos_Serial_Parallel_Range.hpp" #include "initializers.hpp" #include "specfem_setup.hpp" @@ -26,6 +28,106 @@ template struct Power : AnalyticalFunctionType { return std::string("Pow(x,") + std::to_string(power) + ")"; } }; + +/** + * @brief Concatenates k functions into an array-valued function. + * + * This acts just like numpy.concatenate for rank-1 arrays or itertools.chain + * for iterables, creating a vector-valued function, with components coming from + * each `AnalyticalFunctions...` parameter. + * + * @tparam AnalyticalFunctions analytial function types to be concatenated. + */ +template +struct Chain : AnalyticalFunctionType { + static_assert( + ((std::is_base_of_v) && ...), + "Chain needs all of its parameters to be of " + "AnalyticalFunctionType!"); + static constexpr int num_components = + ((AnalyticalFunctions::num_components) + ...); + static std::array + evaluate(const type_real &coord) { + std::array arr; + auto it = arr.begin(); + ( + [&]() { + const auto sub = AnalyticalFunctions::evaluate(coord); + std::copy(sub.begin(), sub.end(), it); + it += AnalyticalFunctions::num_components; + }(), + ...); + return arr; + } + + static std::string description() { + return std::string("Chain (") + std::to_string(num_components) + + " components)\n" + + ((specfem::test_fixture::impl::description::get( + 2) + + "\n") + + ...); + } + static std::string name() { + return std::string("Chain(") + + ((specfem::test_fixture::impl::name::get() + + ",") + + ...) + + ")"; + } +}; + +/** + * @brief Sums k analytical functions. + * + * Represents the pointwise-sum of two or more functions. For example, + * `Sum,AnalyticalFunctionType::Power<0>>` + * represents the function mapping `x` to `x^1 + x^0 = x + 1`. + * + * Each function parameter must have the same dimensions. Unlike `numpy`, we + * have not written any support for broadcasting. + * + * @tparam AnalyticalFunctions analytical function types to be added together. + */ +template struct Sum : AnalyticalFunctionType { + static_assert( + ((std::is_base_of_v) && ...), + "Sum needs all of its parameters to be of " + "AnalyticalFunctionType!"); + static constexpr int num_components = + std::min((AnalyticalFunctions::num_components)...); + + static_assert( + ((AnalyticalFunctions::num_components == num_components) && ...), + "Sum needs all of its parameters to have the same number of components!"); + static std::array + evaluate(const type_real &coord) { + std::array arr{ 0 }; + ( + [&]() { + const auto sub = AnalyticalFunctions::evaluate(coord); + for (int icomp = 0; icomp < num_components; ++icomp) { + arr[icomp] += sub[icomp]; + } + }(), + ...); + return arr; + } + + static std::string description() { + return std::string("Sum (") + std::to_string(num_components) + + " components)\n" + + ((specfem::test_fixture::impl::description::get( + 2) + + "\n") + + ...); + } + static std::string name() { + return std::string("Sum(") + + ((AnalyticalFunctions::description() + ",") + ...) + ")"; + } +}; + } // namespace AnalyticalFunctionType } // namespace specfem::test_fixture diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/initializers.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/initializers.hpp index 2cb2c589f..d8f0e5aa3 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/initializers.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/initializers.hpp @@ -111,4 +111,14 @@ namespace AnalyticalFunctionType { struct AnalyticalFunctionType {}; } // namespace AnalyticalFunctionType +/** + * @brief Provides a struct that can be passed as a LagrangeDerivative argument + * to kernels. + * + * @tparam QuadraturePointsX quadrature on the xi coordinates + * @tparam QuadraturePointsZ quadrature on the gamma coordinates + */ +template +struct LagrangeDerivative2D; + } // namespace specfem::test_fixture diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp new file mode 100644 index 000000000..94eac3b41 --- /dev/null +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp @@ -0,0 +1,49 @@ +#include "../impl/descriptions.hpp" +#include "initializers.hpp" +#include "specfem_setup.hpp" + +namespace specfem::test_fixture { + +template +struct LagrangeDerivative2D { + + static constexpr int ngllx = QuadraturePointsX::nquad; + static constexpr int ngllz = QuadraturePointsZ::nquad; + + using QuadratureRuleX = QuadratureRule; + using QuadratureRuleZ = QuadratureRule; + +private: + using memory_space = Kokkos::DefaultExecutionSpace::memory_space; + +public: + Kokkos::View xi; + Kokkos::View gamma; + +public: + KOKKOS_FUNCTION LagrangeDerivative2D() = default; + KOKKOS_FUNCTION LagrangeDerivative2D(const std::string &name) + : xi(name + "(xi)"), gamma(name + "(gamma)") { + + auto h_xi = Kokkos::create_mirror_view(xi); + auto h_ga = Kokkos::create_mirror_view(gamma); + + for (int ipoly = 0; ipoly < ngllx; ipoly++) { + for (int isample = 0; isample < ngllx; isample++) { + h_xi(isample, ipoly) = QuadratureRuleX::evaluate_lagrange_derivative( + ipoly, QuadraturePointsX::quadrature_points[isample]); + } + } + + for (int ipoly = 0; ipoly < ngllz; ipoly++) { + for (int isample = 0; isample < ngllz; isample++) { + h_ga(isample, ipoly) = QuadratureRuleZ::evaluate_lagrange_derivative( + ipoly, QuadraturePointsZ::quadrature_points[isample]); + } + } + + Kokkos::deep_copy(xi, h_xi); + Kokkos::deep_copy(gamma, h_ga); + } +}; +} // namespace specfem::test_fixture diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp index e37c14f3d..21eb9e077 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp @@ -1,11 +1,12 @@ #pragma once -#include "../impl/descriptions.hpp" +#include "Kokkos_Core.hpp" #include "initializers.hpp" #include "specfem_setup.hpp" #include #include +#include namespace specfem::test_fixture { @@ -13,29 +14,130 @@ template struct QuadratureRule { static_assert(std::is_base_of_v, "QuadratureRule template parameter expects QuadraturePoints!"); - using QuadraturePoints = QuadraturePointsType; + static constexpr int nquad = QuadraturePointsType::nquad; static constexpr std::array quadrature_points = QuadraturePointsType::quadrature_points; - static constexpr type_real evaluate_lagrange_polynomial(const int &iquad, - const type_real &x) { + /** + * @brief Evaluates the `iquad`th Lagrange basis polynomial at the point `x`. + * + * @param iquad - the index of the polynomial (will evaluate to 1 at this knot + * and zero everywhere else) + * @param x - coordinate to evaluate at + * @return double the value L[iquad] (x). + */ + static constexpr double evaluate_lagrange_polynomial(const int &iquad, + const double &x) { + + // TODO: should we switch this to use the polynomial_coefficients? double val = 1; - for (int i = 0; i < nquad; i++) { + for (int i = 0; i < nquad; ++i) { if (i != iquad) { val *= (x - quadrature_points[i]) / (quadrature_points[iquad] - quadrature_points[i]); } } - return (type_real)val; + return val; + } + + /** + * @brief Provides the corresponding lagrange interpolation polynomial over + * the basis {x^k}, that is, $L_{iquad} = \sum_k + * computelagrangepolynomialcoefficients(iquad)[k] * x^k$ + * + * @param iquad index of the lagrange interpolating polynomial to find + * @return std::array the array of coefficients + */ + static constexpr std::array + compute_lagrange_polynomial_coefficients(const int &iquad) { + std::array coeffs{ 0 }; + coeffs[0] = 1; + + for (int i = 0; i < nquad; ++i) { + if (i != iquad) { + // coeffs *= (x - quadrature_points[i])/(quadrature_points[iquad] - + // quadrature_points[i]) + double factor = 1 / (quadrature_points[iquad] - quadrature_points[i]); + coeffs[nquad - 1] = coeffs[nquad - 2] * factor; + for (int j = nquad - 2; j > 0; --j) { + coeffs[j] = + (coeffs[j - 1] - coeffs[j] * quadrature_points[i]) * factor; + } + coeffs[0] *= -quadrature_points[i] * factor; + } + } + return coeffs; + } + + static constexpr std::array, nquad> + _compute_all_lagrange_polynomial_coefficients() { + std::array, nquad> coeffs{}; + for (int i = 0; i < nquad; ++i) { + coeffs[i] = compute_lagrange_polynomial_coefficients(i); + } + return coeffs; + } + + static constexpr std::array, nquad> + polynomial_coefficients = _compute_all_lagrange_polynomial_coefficients(); + + static constexpr double + compute_lagrange_quadrature_weight(const int &iquad, + const double &integral_start = -1, + const double &integral_end = 1) { + std::array L = polynomial_coefficients[iquad]; + + // evaluate integral at integral_end and integral_start + double result = 0; + + double start_pow = integral_start; + double end_pow = integral_end; + for (int i = 0; i < nquad; ++i) { + // add L_integral[k] * integral_end^k - L_integral[k] * integral_start^k + result += (L[i] / (i + 1)) * (end_pow - start_pow); + start_pow *= integral_start; + end_pow *= integral_end; + } + + return result; } - static std::string description(const int &indent = 0) { - return specfem::test_fixture::impl::description::get( - indent); + static constexpr std::array, nquad> + _compute_all_lagrange_polynomial_derivative_coefficients() { + std::array, nquad> coeffs{}; + for (int ipoly = 0; ipoly < nquad; ++ipoly) { + const auto &L = polynomial_coefficients[ipoly]; + for (int ideg = 1; ideg < nquad; ++ideg) { + coeffs[ipoly][ideg - 1] = ideg * L[ideg]; + } + } + return coeffs; } - static std::string quadrature_name() { - return specfem::test_fixture::impl::name::get(); + static constexpr std::array, nquad> + derivative_coefficients = + _compute_all_lagrange_polynomial_derivative_coefficients(); + + /** + * @brief Provides L'(x) for the `iquad`th Lagrange polynomial. + * + * Evaluates the derivative of the `iquad`th Lagrange basis polynomial at the + * point `x`. + * + * @param iquad - the index of the polynomial (L will evaluate to 1 at this + * knot and zero everywhere else) + * @param x - coordinate to evaluate at + * @return double the value L'[iquad] (x). + */ + static constexpr double evaluate_lagrange_derivative(const int &iquad, + const type_real &x) { + double result = 0; + double xpow = 1; + for (int ideg = 0; ideg < nquad - 1; ++ideg) { + result += xpow * derivative_coefficients[iquad][ideg]; + xpow *= x; + } + return result; } }; namespace QuadraturePoints { @@ -44,20 +146,12 @@ struct GLL1 : QuadraturePoints { static constexpr int nquad = 2; static constexpr std::array quadrature_points = { -1, 1 }; - static std::string name() { return "GLL1"; } - static std::string description() { - return ("2-point GLL quadrature (exactness to x^1)\n" - " points = [-1, 1]"); - } + static std::string description() { return "GLL1 (-1, 1)"; } }; struct GLL2 : QuadraturePoints { static constexpr int nquad = 3; static constexpr std::array quadrature_points = { -1, 0, 1 }; - static std::string name() { return "GLL2"; } - static std::string description() { - return ("3-point GLL quadrature (exactness to x^3)\n" - " points = [-1, 0, 1]"); - } + static std::string description() { return "GLL2 (-1, 0, 1)"; } }; struct Asymm5Point : QuadraturePoints { @@ -65,22 +159,18 @@ struct Asymm5Point : QuadraturePoints { static constexpr std::array quadrature_points = { -1, -0.8, -0.5, 0.2, 0.7 }; - static std::string name() { return "Asymm5"; } static std::string description() { - return ("5 point asymmetric quadrature rule (low exactness interpolating " - "quadrature for testing)\n" - " points = [-1, -0.8, -0.5, 0.2, 0.7]"); + return "5 point asymmetric (low exactness interpolating quadrature for " + "testing)"; } }; struct Asymm4Point : QuadraturePoints { static constexpr int nquad = 4; static constexpr std::array quadrature_points = { -0.3, 0, 0.4, 0.6 }; - static std::string name() { return "Asymm4"; } static std::string description() { - return ("4 point asymmetric quadrature rule (low exactness interpolating " - "quadrature for testing)\n" - " points = [-0.3, 0, 0.4, -0.6]"); + return "4 point asymmetric (low exactness interpolating quadrature for " + "testing)"; } }; From 8b2e8e1a188799b99774224a6dff685de9c910b4 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Mon, 19 Jan 2026 18:21:11 -0500 Subject: [PATCH 03/15] completed test -- passes --- .../integrate/coupling_integral1d_dnshape.hpp | 41 +++-- .../dim2/coupling_integral/dshape.cpp | 169 ++++++++++++++---- .../nonconforming_interface/quadrature.hpp | 20 +++ 3 files changed, 178 insertions(+), 52 deletions(-) diff --git a/include/algorithms/integrate/coupling_integral1d_dnshape.hpp b/include/algorithms/integrate/coupling_integral1d_dnshape.hpp index 00cd1f1de..bef1dc8e7 100644 --- a/include/algorithms/integrate/coupling_integral1d_dnshape.hpp +++ b/include/algorithms/integrate/coupling_integral1d_dnshape.hpp @@ -109,9 +109,11 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( // sample derivative of interior_index shape function, in normal // (covector) direction, at edge coordinate. const type_real dlagrange_dn = - ipoint_n_is_ix - ? lagrange_derivative.xi(ipoint_n, self_index.ix) - : lagrange_derivative.gamma(ipoint_n, self_index.iz); + -(ipoint_n_is_ix ? lagrange_derivative.xi(0, ipoint_n) + : lagrange_derivative.gamma(0, ipoint_n)); + // if the edge is at igll = 0, then the outward normal is the negative + // direction. Instead of using self_index.iz and ngll_n - ipoint_n, we + // can use symmetry to assume the edge is at igll == 0 // we may want to refactor this at some point. For now, this should be // fine to update the index point. @@ -145,27 +147,30 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) #pragma unroll #endif + // dshape_dn, local-normal derivative (derivative of + // normal-direction L, which is normally kronecker delta + // indicating on edge) + + // first, in the local-normal component + type_real dshape_dn = + intersection_normal_contravariant_edgelocal(iedge, iquad, 0) * + (dlagrange_dn * transfer_function_self(iquad)); + // dshape_dn, local-tangential derivative (differentiate + // transfer_function_self instead of normal-direction L) + if (ipoint_n == 0) { + // the local-tangential is constant zero if we are not on the + // edge. + dshape_dn += + intersection_normal_contravariant_edgelocal(iedge, iquad, 1) * + transfer_function_self_derivative(iedge, ipoint_s, iquad); + } for (int icomp = 0; icomp < ncomp; icomp++) { - // dshape_dn, local-normal derivative (derivative of - // normal-direction L, which is normally kronecker delta - // indicating on edge) - type_real dshape_dn = - intersection_normal_contravariant_edgelocal(iedge, iquad, 0) * - (dlagrange_dn * transfer_function_self(iquad)); - // dshape_dn, local-tangential derivative (differentiate - // transfer_function_self instead of normal-direction L) - if (ipoint_n == 0) { - dshape_dn += - intersection_normal_contravariant_edgelocal(iedge, iquad, - 1) * - transfer_function_self_derivative(iedge, ipoint_s, iquad); - } result(icomp) += intersection_field(iedge, iquad, icomp) * intersection_factor(iedge, iquad) * dshape_dn; } } - callback(self_index, result); + callback(interior_index, result); } }); } diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index a8ba63b09..86abc930d 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -1,14 +1,14 @@ +#include "Kokkos_Core.hpp" +#include "Kokkos_Core_fwd.hpp" #include "algorithms/integrate/coupling_integral1d_dnshape.hpp" +#include "enumerations/dim2/mesh_entities.hpp" #include "specfem/chunk_edge.hpp" #include "utilities/include/fixture/impl/accessors.hpp" #include "utilities/include/fixture/nonconforming_interface.hpp" #include "SPECFEM_Environment.hpp" -#include "utilities/include/fixture/nonconforming_interface/analytical_function.hpp" -#include "utilities/include/fixture/nonconforming_interface/intersection_function.hpp" #include "utilities/include/fixture/nonconforming_interface/quadrature.hpp" -#include "utilities/include/fixture/nonconforming_interface/transfer_function.hpp" #include template @@ -57,8 +57,8 @@ class ChunkEdgeIterator : public specfem::execution::TeamThreadRangePolicy< return index_type( specfem::point::edge_index( - 0 /*ispec*/, global_offset + iedge /*iedge*/, ipoint, iz, ix, - edge_type), + global_offset + iedge /*ispec*/, global_offset + iedge /*iedge*/, + ipoint, iz, ix, edge_type), iedge, i); } @@ -173,19 +173,22 @@ TEST(CouplingIntegral, SimpleDShapeTest) { using QuadX = specfem::test_fixture::QuadraturePoints::GLL2; using QuadZ = specfem::test_fixture::QuadraturePoints::GLL2; - using QuadIntersection = specfem::test_fixture::QuadraturePoints::Asymm5Point; + using QuadIntersection = specfem::test_fixture::QuadraturePoints::GL6; constexpr int num_edges = 1; const type_real intersection_min = -1; const type_real intersection_max = 1; const auto side = specfem::mesh_entity::dim2::type::top; - // we are integrating d(shape-function)/dn * F, where - // :: F(x) = x^2 - // :: n = [1, x] in the contravariant local coordinate basis (see - // coupling_integral1d_dnshape.hpp) - // :: shape functions on GLL2 x GLL1 elements + // we are integrating d(shape_function)/dn * F, where + // :: F(x) = x + // :: n = [1, x] in the contravariant local edge basis (see + // coupling_integral1d_dnshape.hpp), or [x, 1] in typical local + // coordinates + // :: shape functions on GLL2 x GLL1 elements (shape_function = L_xi(x) * + // L_gamma(z)) // :: integral on [-1,1] (intersection_) // :: edge has same coordinates, and goes between [-1,1] in xi + // :: solution: int(x * (x * L_xi'(x) * L_gamma(z) + L_xi(x) * L_gamma'(z))) // structs to use @@ -195,7 +198,7 @@ TEST(CouplingIntegral, SimpleDShapeTest) { using IntersectionFunctionInitializer = specfem::test_fixture:: IntersectionFunctionInitializer2D::FromAnalyticalFunction< - specfem::test_fixture::AnalyticalFunctionType::Power<2>, + specfem::test_fixture::AnalyticalFunctionType::Power<1>, QuadIntersection>; const auto function_view = specfem::test_fixture::IntersectionFunction2D( IntersectionFunctionInitializer()) @@ -213,9 +216,9 @@ TEST(CouplingIntegral, SimpleDShapeTest) { for (int iedge = 0; iedge < num_edges; iedge++) { for (int iquad = 0; iquad < QuadIntersection::nquad; iquad++) { h_intersection_factor(iedge, iquad) = - specfem::test_fixture::QuadratureRule< - QuadIntersection>::compute_lagrange_quadrature_weight(iquad, -1, - 1); + specfem::test_fixture::QuadratureRule:: + compute_lagrange_quadrature_weight(iquad, intersection_min, + intersection_max); } } intersection_factor.sync_to_device(h_intersection_factor); @@ -236,11 +239,12 @@ TEST(CouplingIntegral, SimpleDShapeTest) { .get_view(); // no data access help for transfer_function_self derivative, but the type - // only needs to support deriv(iedge, iquad, icomp). + // only needs to support deriv(iedge, iedge, iquad). using TransferFunctionDerivativeType = specfem::datatype::VectorChunkEdgeViewType< - type_real, dimension_tag, num_edges, QuadIntersection::nquad, 2, - false, memory_space, Kokkos::MemoryTraits<> >; + type_real, dimension_tag, num_edges, QuadX::nquad, + QuadIntersection::nquad, false, memory_space, + Kokkos::MemoryTraits<> >; TransferFunctionDerivativeType transfer_function_self_derivative( "dshape::transfer_function_self_derivative"); // init later with @@ -320,32 +324,129 @@ TEST(CouplingIntegral, SimpleDShapeTest) { const int num_chunks = num_edges / chunk_size + ((num_edges % chunk_size != 0) ? 1 : 0); + Kokkos::View > + computed_integrals("dshape::computed_integrals", num_edges); + + // Kokkos::View > + // expected_solutions("dshape::expected_solutions", num_edges); + // const auto h_expected_solutions = + // Kokkos::create_mirror_view(expected_solutions); + Kokkos::View + h_expected_solutions("dshape::h_expected_solutions", num_edges); + + // ================================================================= + // compute expected solutions + + using IntersectionQuadrature = + specfem::test_fixture::QuadratureRule; + std::array quadrature_weights; + for (int iquad = 0; iquad < nquad_intersection; iquad++) { + quadrature_weights[iquad] = + IntersectionQuadrature::compute_lagrange_quadrature_weight( + iquad, intersection_min, intersection_max); + } + + // [!] modify for when num_edges > 1 in later test + for (int iz = 0; iz < ngllz; iz++) { + for (int ix = 0; ix < ngllx; ix++) { + // use the intersection quadrature rule + // solution: int(x * (x * L_xi'(x) * L_gamma(z) + L_xi(x) * L_gamma'(z))) + double integral = 0; + for (int iquad = 0; iquad < nquad_intersection; iquad++) { + const double x = QuadIntersection::quadrature_points[iquad]; + const double z = 1; // since we are at the top + const double dshapedxi = + specfem::test_fixture::QuadratureRule< + QuadX>::evaluate_lagrange_derivative(ix, x) * + specfem::test_fixture::QuadratureRule< + QuadZ>::evaluate_lagrange_polynomial(iz, z); + const double dshapedga = + specfem::test_fixture::QuadratureRule< + QuadX>::evaluate_lagrange_polynomial(ix, x) * + specfem::test_fixture::QuadratureRule< + QuadZ>::evaluate_lagrange_derivative(iz, z); + const double intersection_function = + IntersectionFunctionInitializer::AnalyticalFunctionType::evaluate( + x)[0]; + const auto n_contraedge = IntersectionContraNormalFunction::evaluate(x); + double nxi; + double nga; + switch (side) { + case specfem::mesh_entity::dim2::type::bottom: + nga = -n_contraedge[0]; + nxi = n_contraedge[1]; + break; + case specfem::mesh_entity::dim2::type::top: + nga = n_contraedge[0]; + nxi = n_contraedge[1]; + break; + case specfem::mesh_entity::dim2::type::left: + nxi = -n_contraedge[0]; + nga = n_contraedge[1]; + break; + case specfem::mesh_entity::dim2::type::right: + nxi = n_contraedge[0]; + nga = n_contraedge[1]; + break; + default: + FAIL() << "Poorly posed test. \"side\" is not an edge!."; + } + integral += + quadrature_weights[iquad] * + (intersection_function * (nxi * dshapedxi + nga * dshapedga)); + } + h_expected_solutions(0, iz, ix) = (type_real)integral; + } + } + + // Kokkos::deep_copy(expected_solutions, h_expected_solutions); + // ================================================================= + Kokkos::parallel_for( "SimpleDShapeTest", Kokkos::TeamPolicy<>(num_edges, 1, 1), KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type &team_member) { - const int iedge = team_member.league_rank(); - const FunctionType F(Kokkos::subview( - function_view, - Kokkos::make_pair(iedge * chunk_size, (iedge + 1) * chunk_size), - Kokkos::ALL(), Kokkos::ALL())); + const int chunk_index = team_member.league_rank(); + const FunctionType F( + Kokkos::subview(function_view, + Kokkos::make_pair(chunk_index * chunk_size, + (chunk_index + 1) * chunk_size), + Kokkos::ALL(), Kokkos::ALL())); specfem::algorithms::coupling_integral_dnshape( nonconforming_interfaces, ngllz, ngllx, lagrange_derivative, - ChunkEdgeIndex( - Kokkos::subview(edge_types, - Kokkos::make_pair(iedge * chunk_size, - (iedge + 1) * chunk_size)), - num_edges, ngllz, ngllx, team_member), + ChunkEdgeIndex(Kokkos::subview(edge_types, + Kokkos::make_pair( + chunk_index * chunk_size, + (chunk_index + 1) * chunk_size)), + num_edges, ngllz, ngllx, team_member), F, intersection_factor, intersection_contra_normal, transfer_function_self_derivative, [&](const auto &index, const auto &point) { - // TODO START HERE - // for (int icomp = 0; icomp < ncomp_self; ++icomp) { - // Kokkos::single(Kokkos::PerTeam(team_member), [&]() { - // result_view(index(0), index(1), icomp) = point(icomp); - // }); - // } + computed_integrals(index.ispec, index.iz, index.ix) = point(0); }); }); + + const auto h_computed_integrals = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), computed_integrals); + + for (int iedge = 0; iedge < num_edges; iedge++) { + for (int iz = 0; iz < ngllz; iz++) { + for (int ix = 0; ix < ngllx; ix++) { + + if (!specfem::utilities::is_close( + h_computed_integrals(iedge, iz, ix), + h_expected_solutions(iedge, iz, ix))) { + ADD_FAILURE() << "Integral mismatch for edge " << iedge << "\n" + << " at GLL point (iz = " << iz << ", ix = " << ix + << ")\n" + << " expected: " + << h_expected_solutions(iedge, iz, ix) << "\n" + << " computed: " + << h_computed_integrals(iedge, iz, ix); + } + } + } + } } int main(int argc, char *argv[]) { diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp index 21eb9e077..fc9bc5674 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp @@ -173,6 +173,26 @@ struct Asymm4Point : QuadraturePoints { "testing)"; } }; +struct GL6 : QuadraturePoints { + static constexpr int nquad = 7; + + // computed from scipy.special.roots_legendre(7) + // (https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.roots_legendre.html) + // print(*(f"{x:.40}" for x in scipy.special.roots_legendre(7)[0]), sep=", ") + static constexpr std::array quadrature_points = { + -0.9491079123427583752459213428664952516556, + -0.7415311855993944600839995473506860435009, + -0.4058451513773971841558818596240598708391, + 0.0, + 0.4058451513773971841558818596240598708391, + 0.7415311855993944600839995473506860435009, + 0.9491079123427583752459213428664952516556 + }; + static std::string name() { return "GL6"; } + static std::string description() { + return "Degree-6 Gauss-Legendre (11 degrees of exactness)"; + } +}; } // namespace QuadraturePoints From 7b84f12d324bb0d18f0e5b1fc13204c66e50d2cb Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Mon, 19 Jan 2026 21:09:04 -0500 Subject: [PATCH 04/15] 1554 compatibility --- tests/unit-tests/CMakeLists.txt | 8 ++++---- .../algorithms/dim2/coupling_integral/dshape.cpp | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/unit-tests/CMakeLists.txt b/tests/unit-tests/CMakeLists.txt index da86cd740..7b6f2f961 100644 --- a/tests/unit-tests/CMakeLists.txt +++ b/tests/unit-tests/CMakeLists.txt @@ -632,16 +632,16 @@ add_executable( target_link_libraries( coupling_integral_tests - algorithms + specfem::algorithms enumerations - mesh + specfem::mesh assembly gtest_main Kokkos::kokkos specfem_environment - utilities + specfem::utilities - quadrature + specfem::quadrature yaml-cpp -lpthread -lm ) diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index 86abc930d..da693154e 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -1,7 +1,7 @@ #include "Kokkos_Core.hpp" #include "Kokkos_Core_fwd.hpp" -#include "algorithms/integrate/coupling_integral1d_dnshape.hpp" #include "enumerations/dim2/mesh_entities.hpp" +#include "specfem/algorithms/coupling_integral1d_dnshape.hpp" #include "specfem/chunk_edge.hpp" #include "utilities/include/fixture/impl/accessors.hpp" From 7f6ea3de5735265226a441d3079b7962f22d1340 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Mon, 19 Jan 2026 21:11:28 -0500 Subject: [PATCH 05/15] compilation issues on dshape --- .../unit-tests/algorithms/dim2/coupling_integral/dshape.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index da693154e..4332e8405 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -1,7 +1,6 @@ -#include "Kokkos_Core.hpp" -#include "Kokkos_Core_fwd.hpp" -#include "enumerations/dim2/mesh_entities.hpp" #include "specfem/algorithms/coupling_integral1d_dnshape.hpp" + +#include "enumerations/dim2/mesh_entities.hpp" #include "specfem/chunk_edge.hpp" #include "utilities/include/fixture/impl/accessors.hpp" From 52c5a71afd7124630822cee91929f0e350cf4339 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Mon, 19 Jan 2026 23:44:07 -0500 Subject: [PATCH 06/15] purged snuck internal kokkos import --- .../fixture/nonconforming_interface/analytical_function.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp index 837ed2a22..15be5fcd3 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/analytical_function.hpp @@ -1,7 +1,6 @@ #pragma once #include "../impl/descriptions.hpp" -#include "Serial/Kokkos_Serial_Parallel_Range.hpp" #include "initializers.hpp" #include "specfem_setup.hpp" From 611e7cb486580bb6615d296758e95ae187a22fc1 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 20 Jan 2026 01:20:29 -0500 Subject: [PATCH 07/15] fixed unroll pragma location --- core/specfem/algorithms/coupling_integral1d_dnshape.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/core/specfem/algorithms/coupling_integral1d_dnshape.hpp b/core/specfem/algorithms/coupling_integral1d_dnshape.hpp index 39f7bde18..d5c1fbd5f 100644 --- a/core/specfem/algorithms/coupling_integral1d_dnshape.hpp +++ b/core/specfem/algorithms/coupling_integral1d_dnshape.hpp @@ -143,9 +143,6 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( #endif for (int iquad = 0; iquad < nquad_intersection; iquad++) { -#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) -#pragma unroll -#endif // dshape_dn, local-normal derivative (derivative of // normal-direction L, which is normally kronecker delta // indicating on edge) @@ -163,6 +160,9 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( intersection_normal_contravariant_edgelocal(iedge, iquad, 1) * transfer_function_self_derivative(iedge, ipoint_s, iquad); } +#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) +#pragma unroll +#endif for (int icomp = 0; icomp < ncomp; icomp++) { result(icomp) += intersection_field(iedge, iquad, icomp) * intersection_factor(iedge, iquad) * dshape_dn; From 1ccd29988315bb33f24ed1021fc202075c01c2a7 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 20 Jan 2026 15:44:10 -0500 Subject: [PATCH 08/15] fixed gpu compile and bug with accessor axis ordering --- .../algorithms/dim2/coupling_integral/dshape.cpp | 4 +++- .../utilities/include/fixture/impl/accessors.hpp | 15 ++++++++++++++- .../lagrange_derivative.hpp | 2 +- .../nonconforming_interface/quadrature.hpp | 2 +- 4 files changed, 19 insertions(+), 4 deletions(-) diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index 4332e8405..495a9e172 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -160,7 +160,7 @@ class nonconforming_interfaces_patch }; // temporary test for purposes of uncombined coupling_integral -TEST(CouplingIntegral, SimpleDShapeTest) { +void execute_simple_dshape_test() { constexpr auto dimension_tag = specfem::dimension::type::dim2; constexpr auto interface_tag = specfem::interface::interface_tag::acoustic_elastic; @@ -448,6 +448,8 @@ TEST(CouplingIntegral, SimpleDShapeTest) { } } +TEST(CouplingIntegral, SimpleDShapeTest) { execute_simple_dshape_test(); } + int main(int argc, char *argv[]) { ::testing::InitGoogleTest(&argc, argv); ::testing::AddGlobalTestEnvironment(new SPECFEMEnvironment); diff --git a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp index 3959afce7..6c9fca304 100644 --- a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp +++ b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp @@ -32,17 +32,30 @@ struct NonconformingAccessorPatch2D specfem::connections::type::nonconforming; /// View type for storing intersection scaling factors private: + // ====================================================================================== + /** + * @brief helper to unpack [Axes][...] + * + * ViewInternalType::type = T[Axis1]...[Axisk] + */ template struct ViewInternalType; + // stitch in axis recursively. template struct ViewInternalType { + + // I would have thought it would be + // ViewInternalType::type + // but this one gave me the correct type: using type = - typename ViewInternalType::type; + typename ViewInternalType::type[Axis]; }; + // base case (recursive termination) template struct ViewInternalType { using type = LeftType; }; + // ====================================================================================== public: using DataViewType = diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp index 94eac3b41..6b29572be 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp @@ -22,7 +22,7 @@ struct LagrangeDerivative2D { public: KOKKOS_FUNCTION LagrangeDerivative2D() = default; - KOKKOS_FUNCTION LagrangeDerivative2D(const std::string &name) + LagrangeDerivative2D(const std::string &name) : xi(name + "(xi)"), gamma(name + "(gamma)") { auto h_xi = Kokkos::create_mirror_view(xi); diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp index fc9bc5674..8e6beda80 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp @@ -86,7 +86,7 @@ template struct QuadratureRule { compute_lagrange_quadrature_weight(const int &iquad, const double &integral_start = -1, const double &integral_end = 1) { - std::array L = polynomial_coefficients[iquad]; + const std::array &L = polynomial_coefficients[iquad]; // evaluate integral at integral_end and integral_start double result = 0; From 635c8acebf5a4f89eef541f9b300ea85f46dcecb Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 20 Jan 2026 18:44:19 -0500 Subject: [PATCH 09/15] restored removed quadrature name/desc features --- .../include/fixture/impl/accessors.hpp | 2 - .../nonconforming_interface/quadrature.hpp | 41 ++++++++++++++----- 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp index 6c9fca304..9486f5318 100644 --- a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp +++ b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp @@ -1,11 +1,9 @@ #pragma once -#include "Kokkos_Core.hpp" #include "enumerations/coupled_interface.hpp" #include "enumerations/medium.hpp" #include "specfem/data_access/accessor.hpp" #include "specfem/data_access/data_class.hpp" -#include namespace specfem::test_fixture::impl { /** diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp index 8e6beda80..f62841762 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/quadrature.hpp @@ -1,6 +1,6 @@ #pragma once -#include "Kokkos_Core.hpp" +#include "../impl/descriptions.hpp" #include "initializers.hpp" #include "specfem_setup.hpp" @@ -14,7 +14,7 @@ template struct QuadratureRule { static_assert(std::is_base_of_v, "QuadratureRule template parameter expects QuadraturePoints!"); - + using QuadraturePoints = QuadraturePointsType; static constexpr int nquad = QuadraturePointsType::nquad; static constexpr std::array quadrature_points = QuadraturePointsType::quadrature_points; @@ -130,7 +130,7 @@ template struct QuadratureRule { * @return double the value L'[iquad] (x). */ static constexpr double evaluate_lagrange_derivative(const int &iquad, - const type_real &x) { + const double &x) { double result = 0; double xpow = 1; for (int ideg = 0; ideg < nquad - 1; ++ideg) { @@ -139,19 +139,34 @@ template struct QuadratureRule { } return result; } + + static std::string description(const int &indent = 0) { + return specfem::test_fixture::impl::description::get( + indent); + } + static std::string quadrature_name() { + return specfem::test_fixture::impl::name::get(); + } }; namespace QuadraturePoints { struct GLL1 : QuadraturePoints { static constexpr int nquad = 2; static constexpr std::array quadrature_points = { -1, 1 }; - - static std::string description() { return "GLL1 (-1, 1)"; } + static std::string name() { return "GLL1"; } + static std::string description() { + return ("2-point GLL quadrature (exactness to x^1)\n" + " points = [-1, 1]"); + } }; struct GLL2 : QuadraturePoints { static constexpr int nquad = 3; static constexpr std::array quadrature_points = { -1, 0, 1 }; - static std::string description() { return "GLL2 (-1, 0, 1)"; } + static std::string name() { return "GLL2"; } + static std::string description() { + return ("3-point GLL quadrature (exactness to x^3)\n" + " points = [-1, 0, 1]"); + } }; struct Asymm5Point : QuadraturePoints { @@ -159,18 +174,22 @@ struct Asymm5Point : QuadraturePoints { static constexpr std::array quadrature_points = { -1, -0.8, -0.5, 0.2, 0.7 }; + static std::string name() { return "Asymm5"; } static std::string description() { - return "5 point asymmetric (low exactness interpolating quadrature for " - "testing)"; + return ("5 point asymmetric quadrature rule (low exactness interpolating " + "quadrature for testing)\n" + " points = [-1, -0.8, -0.5, 0.2, 0.7]"); } }; struct Asymm4Point : QuadraturePoints { static constexpr int nquad = 4; static constexpr std::array quadrature_points = { -0.3, 0, 0.4, 0.6 }; + static std::string name() { return "Asymm4"; } static std::string description() { - return "4 point asymmetric (low exactness interpolating quadrature for " - "testing)"; + return ("4 point asymmetric quadrature rule (low exactness interpolating " + "quadrature for testing)\n" + " points = [-0.3, 0, 0.4, -0.6]"); } }; struct GL6 : QuadraturePoints { @@ -190,7 +209,7 @@ struct GL6 : QuadraturePoints { }; static std::string name() { return "GL6"; } static std::string description() { - return "Degree-6 Gauss-Legendre (11 degrees of exactness)"; + return "7 Point (degree-6) Gauss-Legendre (11 degrees of exactness)"; } }; From 31a8058fcbd51efdc6b48e7b3a91442603a6c6a3 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Sun, 8 Feb 2026 23:35:32 -0500 Subject: [PATCH 10/15] conflicts post-merge --- .../coupling_integral1d_dnshape.hpp | 15 +++++----- .../dim2/nonconforming_interfaces.hpp | 2 +- tests/unit-tests/CMakeLists.txt | 2 +- .../dim2/coupling_integral/dshape.cpp | 28 +++++++++---------- .../include/fixture/impl/accessors.hpp | 21 +++++++------- 5 files changed, 34 insertions(+), 34 deletions(-) diff --git a/core/specfem/algorithms/coupling_integral1d_dnshape.hpp b/core/specfem/algorithms/coupling_integral1d_dnshape.hpp index d5c1fbd5f..0cc21b94e 100644 --- a/core/specfem/algorithms/coupling_integral1d_dnshape.hpp +++ b/core/specfem/algorithms/coupling_integral1d_dnshape.hpp @@ -1,7 +1,7 @@ #pragma once -#include "enumerations/interface.hpp" #include "specfem/assembly.hpp" +#include "specfem/element_coupling.hpp" #include "specfem/execution.hpp" #include "specfem/point.hpp" #include @@ -52,9 +52,10 @@ namespace specfem::algorithms { * @param callback - callback function to capture integral values * @ingroup AlgorithmsIntegration */ -template KOKKOS_FUNCTION void coupling_integral_dnshape( const specfem::assembly::nonconforming_interfaces @@ -70,7 +71,7 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( , const CallableType &callback) { - constexpr auto self_medium_tag = specfem::interface::attributes< + constexpr auto self_medium_tag = specfem::element_coupling::attributes< dimension_tag, IntersectionFactor::interface_tag>::self_medium(); using PointIndexType = @@ -102,8 +103,8 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( const int ngll_n = ipoint_n_is_ix ? ngllx : ngllz; for (int ipoint_n = 0; ipoint_n < ngll_n; ++ipoint_n) { // iterate backwards from the edge (this gets called for each gllxz) - specfem::point::index interior_index( - self_index.ispec, self_index.iz, self_index.ix); + specfem::point::index + interior_index(self_index.ispec, self_index.iz, self_index.ix); // sample derivative of interior_index shape function, in normal // (covector) direction, at edge coordinate. diff --git a/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp b/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp index 605890a5a..2eac42b50 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim2/nonconforming_interfaces.hpp @@ -20,7 +20,7 @@ class nonconforming_interfaces public: static constexpr auto dimension_tag = specfem::element::dimension_tag::dim2; -private: +protected: template diff --git a/tests/unit-tests/CMakeLists.txt b/tests/unit-tests/CMakeLists.txt index fd22478f8..b2bdf4728 100644 --- a/tests/unit-tests/CMakeLists.txt +++ b/tests/unit-tests/CMakeLists.txt @@ -649,7 +649,7 @@ add_executable( target_link_libraries( coupling_integral_tests specfem::algorithms - enumerations + specfem::element_connections specfem::mesh assembly gtest_main diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index 495a9e172..395c22cc4 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -1,7 +1,7 @@ #include "specfem/algorithms/coupling_integral1d_dnshape.hpp" -#include "enumerations/dim2/mesh_entities.hpp" #include "specfem/chunk_edge.hpp" +#include "specfem/enums.hpp" #include "utilities/include/fixture/impl/accessors.hpp" #include "utilities/include/fixture/nonconforming_interface.hpp" @@ -21,7 +21,7 @@ class ChunkEdgeIterator : public specfem::execution::TeamThreadRangePolicy< specfem::execution::TeamThreadRangePolicy; using execution_space = typename base_type::execution_space; using index_type = - specfem::execution::EdgePointIndex; @@ -55,7 +55,7 @@ class ChunkEdgeIterator : public specfem::execution::TeamThreadRangePolicy< const int ix = is_leftright ? inorm : ipoint; return index_type( - specfem::point::edge_index( + specfem::point::edge_index( global_offset + iedge /*ispec*/, global_offset + iedge /*iedge*/, ipoint, iz, ix, edge_type), iedge, i); @@ -74,7 +74,7 @@ class ChunkEdgeIterator : public specfem::execution::TeamThreadRangePolicy< template class ChunkEdgeIndex { public: static constexpr auto accessor_type = - specfem::data_access::AccessorType::chunk_edge; + specfem::datatype::AccessorType::chunk_edge; using KokkosIndexType = Kokkos::TeamPolicy<>::member_type; using iterator_type = ChunkEdgeIterator; @@ -126,7 +126,7 @@ template class ChunkEdgeIndex { */ class nonconforming_interfaces_patch : public specfem::assembly::nonconforming_interfaces< - specfem::dimension::type::dim2> { + specfem::element::dimension_tag::dim2> { int ngllz; int ngllx; int nquad_intersection; @@ -136,9 +136,9 @@ class nonconforming_interfaces_patch const int &nquad_intersection) : ngllz(ngllz), ngllx(ngllx), nquad_intersection(nquad_intersection) {}; - template + specfem::element_connections::type ConnectionTag> void reinit_container(const int &num_edges) { FOR_EACH_IN_PRODUCT( @@ -161,9 +161,9 @@ class nonconforming_interfaces_patch // temporary test for purposes of uncombined coupling_integral void execute_simple_dshape_test() { - constexpr auto dimension_tag = specfem::dimension::type::dim2; + constexpr auto dimension_tag = specfem::element::dimension_tag::dim2; constexpr auto interface_tag = - specfem::interface::interface_tag::acoustic_elastic; + specfem::element_coupling::interface_tag::acoustic_elastic; constexpr auto boundary_tag = specfem::element::boundary_tag::none; using memory_space = Kokkos::DefaultExecutionSpace::memory_space; @@ -258,21 +258,21 @@ void execute_simple_dshape_test() { constexpr int chunk_size = 1; constexpr auto medium_self = - specfem::interface::attributes::self_medium(); + specfem::element_coupling::attributes::self_medium(); constexpr auto ncomp_self = specfem::element::attributes::components; nonconforming_interfaces_patch nonconforming_interfaces(ngllz, ngllx, nquad_intersection); nonconforming_interfaces.template reinit_container< - interface_tag, boundary_tag, specfem::connections::type::nonconforming>( - num_edges); + interface_tag, boundary_tag, + specfem::element_connections::type::nonconforming>(num_edges); const auto &interface_container = nonconforming_interfaces.template get_interface_container< interface_tag, boundary_tag, - specfem::connections::type::nonconforming>(); + specfem::element_connections::type::nonconforming>(); // ================================================================= // populate this nonconforming interface container and diff --git a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp index 9486f5318..17bce7f7d 100644 --- a/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp +++ b/tests/unit-tests/utilities/include/fixture/impl/accessors.hpp @@ -1,9 +1,8 @@ #pragma once -#include "enumerations/coupled_interface.hpp" -#include "enumerations/medium.hpp" #include "specfem/data_access/accessor.hpp" #include "specfem/data_access/data_class.hpp" +#include "specfem/element_coupling.hpp" namespace specfem::test_fixture::impl { /** @@ -15,19 +14,19 @@ namespace specfem::test_fixture::impl { * @tparam DataClassType * @tparam Axes The size of the view along each axis. */ -template struct NonconformingAccessorPatch2D : specfem::data_access::Accessor< - specfem::data_access::AccessorType::chunk_edge, DataClassType, - specfem::dimension::type::dim2, false /* UseSIMD */> { + specfem::datatype::AccessorType::chunk_edge, DataClassType, + specfem::element::dimension_tag::dim2, false /* UseSIMD */> { public: - static constexpr auto dimension_tag = specfem::dimension::type::dim2; + static constexpr auto dimension_tag = specfem::element::dimension_tag::dim2; static constexpr auto interface_tag = InterfaceTag; static constexpr auto boundary_tag = BoundaryTag; static constexpr auto connection_tag = - specfem::connections::type::nonconforming; + specfem::element_connections::type::nonconforming; /// View type for storing intersection scaling factors private: // ====================================================================================== @@ -90,7 +89,7 @@ struct NonconformingAccessorPatch2D } }; -template struct NonconformingTransferFunctionSelfPatch @@ -107,7 +106,7 @@ struct NonconformingTransferFunctionSelfPatch static constexpr int n_quad_element = NQuadElement; static constexpr int n_quad_intersection = NQuadIntersection; }; -template struct NonconformingTransferFunctionCoupledPatch @@ -124,7 +123,7 @@ struct NonconformingTransferFunctionCoupledPatch static constexpr int n_quad_element = NQuadElement; static constexpr int n_quad_intersection = NQuadIntersection; }; -template struct NonconformingIntersectionNormalPatch @@ -139,7 +138,7 @@ struct NonconformingIntersectionNormalPatch static constexpr int chunk_size = NumberElements; static constexpr int n_quad_intersection = NQuadIntersection; }; -template struct NonconformingIntersectionFactorPatch From 27ec141256bbf128a3ba03ef0f1adb0ac66be422 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Mon, 9 Feb 2026 00:28:49 -0500 Subject: [PATCH 11/15] shape function deriv precompute kernel start --- .../shape_function_normal_derivative.hpp | 139 ++++++++++++++++++ .../dim2/coupling_integral/dshape.cpp | 7 + 2 files changed, 146 insertions(+) create mode 100644 core/specfem/algorithms/shape_function_normal_derivative.hpp diff --git a/core/specfem/algorithms/shape_function_normal_derivative.hpp b/core/specfem/algorithms/shape_function_normal_derivative.hpp new file mode 100644 index 000000000..f44a880cf --- /dev/null +++ b/core/specfem/algorithms/shape_function_normal_derivative.hpp @@ -0,0 +1,139 @@ +#pragma once + +#include "specfem/assembly.hpp" +#include "specfem/assembly/edge_types.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/element/tags.hpp" +#include "specfem/element_coupling.hpp" +#include "specfem/execution.hpp" +#include "specfem/point.hpp" +#include + +namespace specfem::algorithms { + +/** + * @brief Takes a field `intersection_field` on the intersection and computes, + * for each self GLL point, the integral of `intersection_field` times the + * normal derivative of the shape function at that point. `intersection_field` + * should be call-accessible (e.g. Kokkos::View) with shape: + * + * (chunk_size, n_quad_intersection, self::components()) + * + * After handling any other intersection forces, boundary conditions, etc. the + * result can be `atomic_add`ed to the acceleration field. + * + * @tparam dimension_tag dimension of the simulation + * @param nonconforming_interfaces - assembly.nonconforming_interfaces struct + */ +template +Kokkos::View +shape_function_self_normal_derivatives( + const specfem::assembly::edge_types &edge_types, + const specfem::assembly::mesh &mesh, + const specfem::assembly::nonconforming_interfaces + &nonconforming_interfaces) { + using ReturnViewType = + Kokkos::View; + + const auto [self_edges, coupled_edges] = edge_types.get_edges_on_device( + specfem::element_connections::type::nonconforming, interface_tag, + boundary_tag); + const auto &element_grid = mesh.element_grid; + + // TODO: get nquad_intersection from somewhere else + const int ngllz = element_grid.ngllz; + const int ngllx = element_grid.ngllx; + const int nquad_intersection = std::max(ngllz, ngllx); + ReturnViewType normal_derivs("shape_function_self_normal_derivatives", + self_edges.n_edges, ngllz, ngllx, + nquad_intersection); + + // using ParallelConfig = + // specfem::parallel_configuration::default_chunk_config< + // dimension_tag, specfem::datatype::simd, Kokkos::DefaultExecutionSpace>; + // specfem::execution::ChunkedDomainIterator chunk( + // ParallelConfig(), self_edges.element_index, element_grid); + + using parallel_config = + specfem::parallel_configuration::default_chunk_edge_config< + dimension_tag, Kokkos::DefaultExecutionSpace>; + specfem::execution::ChunkedEdgeIterator chunk(parallel_config(), self_edges); + + return normal_derivs; // remove when done + + specfem::execution::for_each_level( + "specfem::compute::compute_stiffness_interaction", chunk, + KOKKOS_LAMBDA( + const typename decltype(chunk)::index_type &chunk_iterator_index) { + const auto &chunk_index = chunk_iterator_index.get_index(); + const auto &team = chunk_index.get_policy_index(); + const int &num_edges = chunk_index.nedges(); + + specfem::execution::for_each_level( + specfem::execution::TeamThreadMDRangeIterator(team, num_edges, + ngllz, ngllx), + [&](const auto &index) { + const int iedge = index(0); + const int iz = index(1); + const int ix = index(2); + + // this index will always have ipoint = 0, but we will not use it + const auto edge_index = + chunk_index.get_iterator()(iedge).get_index(); + + const bool ipoint_n_is_ix = + edge_index.edge_type == + specfem::mesh_entity::dim2::type::left || + edge_index.edge_type == + specfem::mesh_entity::dim2::type::right; + const int ipoint_s = ipoint_n_is_ix ? iz : ix; + const int ipoint_n = ipoint_n_is_ix ? ix : iz; + + // TODO find access to lagrange_derivative and uncomment + const type_real dlagrange_dn = 0; + // = -(ipoint_n_is_ix ? lagrange_derivative.xi(0, ipoint_n) + // : lagrange_derivative.gamma(0, + // ipoint_n)); + + for (int iquad = 0; iquad < nquad_intersection; iquad++) { + + // dshape_dn, local-normal derivative (derivative of + // normal-direction L, which is normally kronecker delta + // indicating on edge) + + // first, in the local-normal component + type_real dshape_dn = 0; + // = intersection_normal_contravariant_edgelocal(iedge, iquad, + // 0) * (dlagrange_dn * transfer_function_self(iquad)); + // dshape_dn, local-tangential derivative (differentiate + // transfer_function_self instead of normal-direction L) + if (ipoint_n == 0) { + // the local-tangential is constant zero if we are not on the + // edge. + // dshape_dn += + // intersection_normal_contravariant_edgelocal(iedge, + // iquad, 1) * transfer_function_self_derivative(iedge, + // ipoint_s, iquad); + } + + normal_derivs(edge_index.iedge, iz, ix, iquad) = dshape_dn; + } + }); + // specfem::execution::for_each_level( + // chunk_index.get_iterator(), + // [&](const typename ChunkIndexType::iterator_type::index_type + // &iterator_index) { + // const auto index = iterator_index.get_index(); + // const auto index_local = iterator_index.get_local_index(); + + // index.ielem + // }); + }); + + return normal_derivs; +} + +} // namespace specfem::algorithms diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index 395c22cc4..c7071a24e 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -1,4 +1,5 @@ #include "specfem/algorithms/coupling_integral1d_dnshape.hpp" +#include "specfem/algorithms/shape_function_normal_derivative.hpp" #include "specfem/chunk_edge.hpp" #include "specfem/enums.hpp" @@ -402,6 +403,12 @@ void execute_simple_dshape_test() { // Kokkos::deep_copy(expected_solutions, h_expected_solutions); // ================================================================= + specfem::assembly::edge_types assembly_edge_types; + specfem::assembly::mesh mesh; + specfem::algorithms::shape_function_self_normal_derivatives( + assembly_edge_types, mesh, nonconforming_interfaces); + Kokkos::parallel_for( "SimpleDShapeTest", Kokkos::TeamPolicy<>(num_edges, 1, 1), KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type &team_member) { From 20f425941a0a5182261650b1c5e35714995ee99f Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Mon, 9 Feb 2026 16:27:20 -0500 Subject: [PATCH 12/15] completed precompute --- .../coupling_integral1d_dnshape.hpp | 147 +++++------------- .../shape_function_normal_derivative.hpp | 94 +++++++---- .../dim2/coupling_integral/dshape.cpp | 69 +++++++- 3 files changed, 167 insertions(+), 143 deletions(-) diff --git a/core/specfem/algorithms/coupling_integral1d_dnshape.hpp b/core/specfem/algorithms/coupling_integral1d_dnshape.hpp index 0cc21b94e..92a40f5c2 100644 --- a/core/specfem/algorithms/coupling_integral1d_dnshape.hpp +++ b/core/specfem/algorithms/coupling_integral1d_dnshape.hpp @@ -26,50 +26,39 @@ namespace specfem::algorithms { * After handling any other intersection forces, boundary conditions, etc. the * result can be `atomic_add`ed to the acceleration field. * - * @tparam dimension_tag dimension of the simulation * @tparam IndexType The chunk_edge iterator type * @tparam IntersectionFieldViewType The type of `intersection_field` - * @tparam ChunkEdgeWeightJacobianType A nonconforming chunk_edge accessor + * @tparam IntersectionFactor A nonconforming chunk_edge accessor * holding `intersection_factor` + * @tparam ShapeFunctionNormalDerivativesType a view holding the normal + * derivatives of the shape function (implementation may change) * @tparam CallableType The callback function, which will be given the point * index and corresponding evaluated integral - * @param nonconforming_interfaces - assembly.nonconforming_interfaces struct - * @param lagrange_derivative - Local derivatives of shape functions on the - * element + * + * @param ngllz - number of quadrature points in z + * @param ngllx - number of quadrature points in x * @param chunk_index - the outer index (chunk_edge) that gets iterated for * points * @param intersection_field - the field to integrate - * @param weight_jacobian - nonconforming chunk_edge accessor holding + * @param intersection_factor - nonconforming chunk_edge accessor holding * `intersection_factor` - * @param intersection_normal_contravariant_edgelocal - the normal vector at - * each intersection point in contravariant local coordinate basis, with the - * first component in the outward normal direction (-/+xi for left/right, - * +/-gamma for top/bottom), and second component in the tangential axis (+gamma - * for left/right, +xi for top/bottom). - * @param transfer_function_self_derivative - (TEMPORARY, TO BE LOADED LATER BY - * THIS FUNCTION) derivative (in edge coordinates, not intersection coordinates - * ) of transfer_function_self, evaluated at intersection quadrature points. + * @param shape_function_normal_derivatives - (TEMPORARY, TO BE LOADED LATER BY + * THIS FUNCTION) normal derivatives of the shape functions at each intersection + * point. * @param callback - callback function to capture integral values * @ingroup AlgorithmsIntegration */ -template +template KOKKOS_FUNCTION void coupling_integral_dnshape( - const specfem::assembly::nonconforming_interfaces - &nonconforming_interfaces, - const int &ngllz, const int &ngllx, - const QuadratureType &lagrange_derivative, const IndexType &chunk_index, + const int &ngllz, const int &ngllx, const IndexType &chunk_index, const IntersectionFieldViewType &intersection_field, const IntersectionFactor &intersection_factor, - const IntersectionJacobianType &intersection_normal_contravariant_edgelocal, - const TransferFunctionDerivativeType - &transfer_function_self_derivative /* This could be loaded just like - transfer_function_self */ - , + const ShapeFunctionNormalDerivativesType &shape_function_normal_derivatives, const CallableType &callback) { + constexpr specfem::element::dimension_tag dimension_tag = + specfem::element::dimension_tag::dim2; constexpr auto self_medium_tag = specfem::element_coupling::attributes< dimension_tag, IntersectionFactor::interface_tag>::self_medium(); @@ -79,99 +68,49 @@ KOKKOS_FUNCTION void coupling_integral_dnshape( using PointFieldType = specfem::point::acceleration; - using SelfTransferFunctionType = specfem::point::transfer_function_self< - IntersectionFactor::n_quad_intersection, dimension_tag, - IntersectionFactor::interface_tag, IntersectionFactor::boundary_tag>; constexpr int ncomp = PointFieldType::components; constexpr int nquad_intersection = IntersectionFactor::n_quad_intersection; + const auto &team = chunk_index.get_policy_index(); + const int &num_edges = chunk_index.nedges(); specfem::execution::for_each_level( - chunk_index.get_iterator(), - [&](const typename IndexType::iterator_type::index_type &index) { - const auto self_index = index.get_index(); - const auto self_index_local = index.get_local_index(); - const bool ipoint_n_is_ix = - self_index.edge_type == specfem::mesh_entity::dim2::type::left || - self_index.edge_type == specfem::mesh_entity::dim2::type::right; - const int &iedge = self_index_local.iedge; - SelfTransferFunctionType transfer_function_self; - specfem::assembly::load_on_device(self_index, nonconforming_interfaces, - transfer_function_self); - - const int &ipoint_s = self_index.ipoint; - const int ngll_n = ipoint_n_is_ix ? ngllx : ngllz; - for (int ipoint_n = 0; ipoint_n < ngll_n; ++ipoint_n) { - // iterate backwards from the edge (this gets called for each gllxz) - specfem::point::index - interior_index(self_index.ispec, self_index.iz, self_index.ix); - - // sample derivative of interior_index shape function, in normal - // (covector) direction, at edge coordinate. - const type_real dlagrange_dn = - -(ipoint_n_is_ix ? lagrange_derivative.xi(0, ipoint_n) - : lagrange_derivative.gamma(0, ipoint_n)); - // if the edge is at igll = 0, then the outward normal is the negative - // direction. Instead of using self_index.iz and ngll_n - ipoint_n, we - // can use symmetry to assume the edge is at igll == 0 - - // we may want to refactor this at some point. For now, this should be - // fine to update the index point. - switch (self_index.edge_type) { - case specfem::mesh_entity::dim2::type::right: - interior_index.ix = (ngllx - 1) - ipoint_n; - break; - case specfem::mesh_entity::dim2::type::top: - interior_index.iz = (ngllz - 1) - ipoint_n; - break; - case specfem::mesh_entity::dim2::type::left: - interior_index.ix = ipoint_n; - break; - case specfem::mesh_entity::dim2::type::bottom: - interior_index.iz = ipoint_n; - break; - } - - PointFieldType result; + specfem::execution::TeamThreadMDRangeIterator(team, num_edges, ngllz, + ngllx), + [&](const auto &index) { + const int iedge = index(0); + const int iz = index(1); + const int ix = index(2); + + PointFieldType result; #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) #pragma unroll #endif - for (int icomp = 0; icomp < ncomp; icomp++) { - result(icomp) = 0; - } + for (int icomp = 0; icomp < ncomp; icomp++) { + result(icomp) = 0; + } + #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) #pragma unroll #endif - for (int iquad = 0; iquad < nquad_intersection; iquad++) { - - // dshape_dn, local-normal derivative (derivative of - // normal-direction L, which is normally kronecker delta - // indicating on edge) + for (int iquad = 0; iquad < nquad_intersection; iquad++) { - // first, in the local-normal component - type_real dshape_dn = - intersection_normal_contravariant_edgelocal(iedge, iquad, 0) * - (dlagrange_dn * transfer_function_self(iquad)); - // dshape_dn, local-tangential derivative (differentiate - // transfer_function_self instead of normal-direction L) - if (ipoint_n == 0) { - // the local-tangential is constant zero if we are not on the - // edge. - dshape_dn += - intersection_normal_contravariant_edgelocal(iedge, iquad, 1) * - transfer_function_self_derivative(iedge, ipoint_s, iquad); - } #if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) #pragma unroll #endif - for (int icomp = 0; icomp < ncomp; icomp++) { - result(icomp) += intersection_field(iedge, iquad, icomp) * - intersection_factor(iedge, iquad) * dshape_dn; - } + for (int icomp = 0; icomp < ncomp; icomp++) { + result(icomp) += + intersection_field(iedge, iquad, icomp) * + intersection_factor(iedge, iquad) * + shape_function_normal_derivatives(iedge, iz, ix, iquad); } - - callback(interior_index, result); } + // this index will always have ipoint = 0, but we will not use it + const auto edge_index = chunk_index.get_iterator()(iedge).get_index(); + + specfem::point::index self_index(edge_index.ispec, iz, + ix); + callback(self_index, result); }); } diff --git a/core/specfem/algorithms/shape_function_normal_derivative.hpp b/core/specfem/algorithms/shape_function_normal_derivative.hpp index f44a880cf..32b685d6d 100644 --- a/core/specfem/algorithms/shape_function_normal_derivative.hpp +++ b/core/specfem/algorithms/shape_function_normal_derivative.hpp @@ -3,6 +3,7 @@ #include "specfem/assembly.hpp" #include "specfem/assembly/edge_types.hpp" #include "specfem/assembly/mesh.hpp" +#include "specfem/chunk_edge/nonconforming_interface.hpp" #include "specfem/element/tags.hpp" #include "specfem/element_coupling.hpp" #include "specfem/execution.hpp" @@ -26,26 +27,31 @@ namespace specfem::algorithms { * @param nonconforming_interfaces - assembly.nonconforming_interfaces struct */ template + specfem::element::boundary_tag boundary_tag, int nquad_element_, + int nquad_intersection_, typename SelfEdgeListView, + typename HPrimeView, typename TransformedIntersectionNormalView> Kokkos::View shape_function_self_normal_derivatives( - const specfem::assembly::edge_types &edge_types, - const specfem::assembly::mesh &mesh, - const specfem::assembly::nonconforming_interfaces - &nonconforming_interfaces) { + const SelfEdgeListView &self_edges, + const specfem::assembly::nonconforming_interfaces< + specfem::element::dimension_tag::dim2> &nonconforming_interfaces, + const int &ngllz, const int &ngllx, const HPrimeView &hprime, + const TransformedIntersectionNormalView + &intersection_normal_contravariant_edgelocal) { + const auto dimension_tag = specfem::element::dimension_tag::dim2; using ReturnViewType = Kokkos::View; - const auto [self_edges, coupled_edges] = edge_types.get_edges_on_device( - specfem::element_connections::type::nonconforming, interface_tag, - boundary_tag); - const auto &element_grid = mesh.element_grid; - // TODO: get nquad_intersection from somewhere else - const int ngllz = element_grid.ngllz; - const int ngllx = element_grid.ngllx; - const int nquad_intersection = std::max(ngllz, ngllx); + const int nquad_intersection = nquad_intersection_; + if (nquad_intersection != nquad_intersection_) { + throw std::runtime_error( + std::string("shape_function_self_normal_derivatives() kernel run with " + "nquad_intersection = ") + + std::to_string(nquad_intersection_) + + ", but assembly got nquad_intersection == " + + std::to_string(nquad_intersection)); + } ReturnViewType normal_derivs("shape_function_self_normal_derivatives", self_edges.n_edges, ngllz, ngllx, nquad_intersection); @@ -62,16 +68,24 @@ shape_function_self_normal_derivatives( dimension_tag, Kokkos::DefaultExecutionSpace>; specfem::execution::ChunkedEdgeIterator chunk(parallel_config(), self_edges); - return normal_derivs; // remove when done + using TransferFunctionSelf = specfem::chunk_edge::transfer_function_self< + dimension_tag, interface_tag, boundary_tag, parallel_config::chunk_size, + nquad_intersection_, nquad_element_>; specfem::execution::for_each_level( - "specfem::compute::compute_stiffness_interaction", chunk, + "specfem::compute::shape_function_normal_derivative", + chunk.set_scratch_size( + 0, Kokkos::PerTeam(TransferFunctionSelf::shmem_size())), KOKKOS_LAMBDA( const typename decltype(chunk)::index_type &chunk_iterator_index) { const auto &chunk_index = chunk_iterator_index.get_index(); const auto &team = chunk_index.get_policy_index(); const int &num_edges = chunk_index.nedges(); + TransferFunctionSelf transfer_function_self(team); + specfem::assembly::load_on_device(chunk_index, nonconforming_interfaces, + transfer_function_self); + specfem::execution::for_each_level( specfem::execution::TeamThreadMDRangeIterator(team, num_edges, ngllz, ngllx), @@ -90,33 +104,51 @@ shape_function_self_normal_derivatives( edge_index.edge_type == specfem::mesh_entity::dim2::type::right; const int ipoint_s = ipoint_n_is_ix ? iz : ix; - const int ipoint_n = ipoint_n_is_ix ? ix : iz; + int ipoint_n = ipoint_n_is_ix ? ix : iz; + const int ngll_n = ipoint_n_is_ix ? ngllx : ngllz; + const int ngll_s = ipoint_n_is_ix ? ngllz : ngllx; + + if (edge_index.edge_type == + specfem::mesh_entity::dim2::type::right || + edge_index.edge_type == + specfem::mesh_entity::dim2::type::top) { + // we want 0 to be on the edge + ipoint_n = ngll_n - 1 - ipoint_n; + } - // TODO find access to lagrange_derivative and uncomment - const type_real dlagrange_dn = 0; - // = -(ipoint_n_is_ix ? lagrange_derivative.xi(0, ipoint_n) - // : lagrange_derivative.gamma(0, - // ipoint_n)); + const type_real dlagrange_dn = -hprime(0, ipoint_n); for (int iquad = 0; iquad < nquad_intersection; iquad++) { // dshape_dn, local-normal derivative (derivative of - // normal-direction L, which is normally kronecker delta - // indicating on edge) + // normal-direction L, which used to be kronecker delta + // indicating on edge. Now has dlagrange_dn != 0 for all + // ipoint_n) // first, in the local-normal component - type_real dshape_dn = 0; - // = intersection_normal_contravariant_edgelocal(iedge, iquad, - // 0) * (dlagrange_dn * transfer_function_self(iquad)); + type_real dshape_dn = + intersection_normal_contravariant_edgelocal( + edge_index.iedge, iquad, 0) * + (dlagrange_dn * + transfer_function_self(iedge, ipoint_s, iquad)); // dshape_dn, local-tangential derivative (differentiate // transfer_function_self instead of normal-direction L) if (ipoint_n == 0) { // the local-tangential is constant zero if we are not on the // edge. - // dshape_dn += - // intersection_normal_contravariant_edgelocal(iedge, - // iquad, 1) * transfer_function_self_derivative(iedge, - // ipoint_s, iquad); + + // recover L' (at point ipoint_s) on intersection quadrature + // points: use transfer function. + type_real transfer_function_derivative = 0; + for (int ipoint_s2 = 0; ipoint_s2 < ngll_s; ipoint_s2++) { + transfer_function_derivative += + hprime(ipoint_s2, ipoint_s) * + transfer_function_self(iedge, ipoint_s2, iquad); + } + + dshape_dn += intersection_normal_contravariant_edgelocal( + edge_index.iedge, iquad, 1) * + transfer_function_derivative; } normal_derivs(edge_index.iedge, iz, ix, iquad) = dshape_dn; diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index c7071a24e..9982c61de 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -315,6 +315,34 @@ void execute_simple_dshape_test() { } Kokkos::deep_copy(edge_types, h_edge_types); + specfem::assembly::EdgeView edgelist( + "dshape::edgelist", num_edges, nquad_element); + { + specfem::assembly::EdgeView h_edgelist( + "dshape::h_edgelist", num_edges, nquad_element); + const auto element = specfem::mesh_entity::element(ngllz, ngllx); + for (int iedge = 0; iedge < num_edges; iedge++) { + h_edgelist.element_index(iedge) = iedge; + h_edgelist.edge_index(iedge) = iedge; + h_edgelist.edge_types(iedge) = h_edge_types(iedge); + + for (int ipoint = 0; ipoint < nquad_element; ipoint++) { + const auto [iz, ix] = + element.map_coordinates(h_edge_types(iedge), ipoint); + h_edgelist.iz(iedge, ipoint) = iz; + h_edgelist.ix(iedge, ipoint) = ix; + } + } + + // specfem::assembly::edge_types::deep_copy(edgelist, + // h_edgelist); + Kokkos::deep_copy(edgelist.element_index, h_edgelist.element_index); + Kokkos::deep_copy(edgelist.edge_index, h_edgelist.edge_index); + Kokkos::deep_copy(edgelist.edge_types, h_edgelist.edge_types); + Kokkos::deep_copy(edgelist.iz, h_edgelist.iz); + Kokkos::deep_copy(edgelist.ix, h_edgelist.ix); + } + // chunk subviews (FunctionType comes from copied test template -- skip the // rest and assume num_edges == 1) using FunctionType = specfem::datatype::VectorChunkEdgeViewType< @@ -403,11 +431,19 @@ void execute_simple_dshape_test() { // Kokkos::deep_copy(expected_solutions, h_expected_solutions); // ================================================================= - specfem::assembly::edge_types assembly_edge_types; - specfem::assembly::mesh mesh; - specfem::algorithms::shape_function_self_normal_derivatives( - assembly_edge_types, mesh, nonconforming_interfaces); + constexpr int nquad_element_ = 3; + constexpr int nquad_intersection_ = 7; + if (nquad_element != nquad_element_) { + throw std::runtime_error("dshape: Wrong kernel for nquad_element!"); + } + if (nquad_intersection != nquad_intersection_) { + throw std::runtime_error("dshape: Wrong kernel for nquad_intersection!"); + } + const auto shape_function_normal_derivatives = + specfem::algorithms::shape_function_self_normal_derivatives< + interface_tag, boundary_tag, nquad_element_, nquad_intersection_>( + edgelist, nonconforming_interfaces, ngllz, ngllx, + lagrange_derivative.xi, intersection_contra_normal); Kokkos::parallel_for( "SimpleDShapeTest", Kokkos::TeamPolicy<>(num_edges, 1, 1), @@ -418,18 +454,35 @@ void execute_simple_dshape_test() { Kokkos::make_pair(chunk_index * chunk_size, (chunk_index + 1) * chunk_size), Kokkos::ALL(), Kokkos::ALL())); + const auto SFND = + Kokkos::subview(shape_function_normal_derivatives, + Kokkos::make_pair(chunk_index * chunk_size, + (chunk_index + 1) * chunk_size), + Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); specfem::algorithms::coupling_integral_dnshape( - nonconforming_interfaces, ngllz, ngllx, lagrange_derivative, + ngllz, ngllx, ChunkEdgeIndex(Kokkos::subview(edge_types, Kokkos::make_pair( chunk_index * chunk_size, (chunk_index + 1) * chunk_size)), num_edges, ngllz, ngllx, team_member), - F, intersection_factor, intersection_contra_normal, - transfer_function_self_derivative, + F, intersection_factor, SFND, [&](const auto &index, const auto &point) { computed_integrals(index.ispec, index.iz, index.ix) = point(0); }); + // specfem::algorithms::coupling_integral_dnshape( + // nonconforming_interfaces, ngllz, ngllx, lagrange_derivative, + // ChunkEdgeIndex(Kokkos::subview(edge_types, + // Kokkos::make_pair( + // chunk_index * chunk_size, + // (chunk_index + 1) * + // chunk_size)), + // num_edges, ngllz, ngllx, team_member), + // F, intersection_factor, intersection_contra_normal, + // transfer_function_self_derivative, + // [&](const auto &index, const auto &point) { + // computed_integrals(index.ispec, index.iz, index.ix) = point(0); + // }); }); const auto h_computed_integrals = Kokkos::create_mirror_view_and_copy( From beda08825b60211cdadcc7916aa644c289c70a50 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 17 Feb 2026 09:21:18 -0500 Subject: [PATCH 13/15] removed iterator patch, moistened edgeview code --- .../shape_function_normal_derivative.hpp | 26 +-- .../assembly/edge_types/dim2/edge_types.cpp | 57 +++--- .../assembly/edge_types/dim2/edge_types.hpp | 10 ++ .../dim2/coupling_integral/dshape.cpp | 166 +++--------------- 4 files changed, 85 insertions(+), 174 deletions(-) diff --git a/core/specfem/algorithms/shape_function_normal_derivative.hpp b/core/specfem/algorithms/shape_function_normal_derivative.hpp index 32b685d6d..9481e753f 100644 --- a/core/specfem/algorithms/shape_function_normal_derivative.hpp +++ b/core/specfem/algorithms/shape_function_normal_derivative.hpp @@ -56,13 +56,6 @@ shape_function_self_normal_derivatives( self_edges.n_edges, ngllz, ngllx, nquad_intersection); - // using ParallelConfig = - // specfem::parallel_configuration::default_chunk_config< - // dimension_tag, specfem::datatype::simd, Kokkos::DefaultExecutionSpace>; - // specfem::execution::ChunkedDomainIterator chunk( - // ParallelConfig(), self_edges.element_index, element_grid); - using parallel_config = specfem::parallel_configuration::default_chunk_edge_config< dimension_tag, Kokkos::DefaultExecutionSpace>; @@ -98,21 +91,34 @@ shape_function_self_normal_derivatives( const auto edge_index = chunk_index.get_iterator()(iedge).get_index(); + /* + * We convert local coordinates / indices ix, iz + * into an edge-local coordinate system: + * + * ================================== ipoint_s = ix for top + * | and bottom, iz for left + * | and right + * | + * ipoint_n = ngll-1 + * + * ipoint_n = 0 represents the edge. + */ + const bool ipoint_n_is_ix = edge_index.edge_type == specfem::mesh_entity::dim2::type::left || edge_index.edge_type == specfem::mesh_entity::dim2::type::right; - const int ipoint_s = ipoint_n_is_ix ? iz : ix; - int ipoint_n = ipoint_n_is_ix ? ix : iz; const int ngll_n = ipoint_n_is_ix ? ngllx : ngllz; const int ngll_s = ipoint_n_is_ix ? ngllz : ngllx; + const int ipoint_s = ipoint_n_is_ix ? iz : ix; + int ipoint_n = ipoint_n_is_ix ? ix : iz; if (edge_index.edge_type == specfem::mesh_entity::dim2::type::right || edge_index.edge_type == specfem::mesh_entity::dim2::type::top) { - // we want 0 to be on the edge + // flip ipoint_n coord so edge is 0 ipoint_n = ngll_n - 1 - ipoint_n; } diff --git a/core/specfem/assembly/edge_types/dim2/edge_types.cpp b/core/specfem/assembly/edge_types/dim2/edge_types.cpp index 0777e96a4..c001ff4a0 100644 --- a/core/specfem/assembly/edge_types/dim2/edge_types.cpp +++ b/core/specfem/assembly/edge_types/dim2/edge_types.cpp @@ -98,28 +98,13 @@ specfem::assembly::edge_types:: "specfem::assembly::interface_types::self_edges", count, ngll); _coupled_edges_ = EdgeViewType( "specfem::assembly::interface_types::coupled_edges", count, ngll); - _h_self_edges_ = edge_types::create_mirror_view(_self_edges_); - _h_coupled_edges_ = edge_types::create_mirror_view(_coupled_edges_); - - for (int iedge = 0; iedge < count; iedge++) { - _h_self_edges_.element_index(iedge) = self_collect[iedge].ispec; - _h_self_edges_.edge_index(iedge) = self_collect[iedge].iedge; - _h_self_edges_.edge_types(iedge) = self_collect[iedge].edge_type; - _h_coupled_edges_.element_index(iedge) = coupled_collect[iedge].ispec; - _h_coupled_edges_.edge_index(iedge) = coupled_collect[iedge].iedge; - _h_coupled_edges_.edge_types(iedge) = - coupled_collect[iedge].edge_type; - for (int ipoint = 0; ipoint < ngll; ipoint++) { - const auto [iz1, ix1] = - element.map_coordinates(self_collect[iedge].edge_type, ipoint); - _h_self_edges_.iz(iedge, ipoint) = iz1; - _h_self_edges_.ix(iedge, ipoint) = ix1; - const auto [iz2, ix2] = element.map_coordinates( - coupled_collect[iedge].edge_type, ipoint); - _h_coupled_edges_.iz(iedge, ipoint) = iz2; - _h_coupled_edges_.ix(iedge, ipoint) = ix2; - } - } + + _h_self_edges_ = edge_view_from_collected_edges( + "specfem::assembly::interface_types::self_edges_host_mirror", + self_collect, element); + _h_coupled_edges_ = edge_view_from_collected_edges( + "specfem::assembly::interface_types::coupled_edges_host_mirror", + coupled_collect, element); edge_types::deep_copy(_self_edges_, _h_self_edges_); edge_types::deep_copy(_coupled_edges_, _h_coupled_edges_); @@ -171,3 +156,31 @@ specfem::assembly::edge_types:: throw std::runtime_error( "Connection type, interface type or boundary type not found"); } + +specfem::assembly::edge_types< + specfem::element::dimension_tag::dim2>::EdgeViewType::HostMirror +specfem::assembly::edge_view_from_collected_edges( + const std::string &label, + const std::vector< + specfem::mesh_entity::edge > + &self_collect, + const specfem::mesh_entity::element + &element) { + const int &ngll = element.ngllx; + const int &count = self_collect.size(); + + specfem::assembly::edge_types:: + EdgeViewType::HostMirror self_edges(label, count, ngll); + for (int iedge = 0; iedge < count; iedge++) { + self_edges.element_index(iedge) = self_collect[iedge].ispec; + self_edges.edge_index(iedge) = self_collect[iedge].iedge; + self_edges.edge_types(iedge) = self_collect[iedge].edge_type; + for (int ipoint = 0; ipoint < ngll; ipoint++) { + const auto [iz1, ix1] = + element.map_coordinates(self_collect[iedge].edge_type, ipoint); + self_edges.iz(iedge, ipoint) = iz1; + self_edges.ix(iedge, ipoint) = ix1; + } + } + return self_edges; +} diff --git a/core/specfem/assembly/edge_types/dim2/edge_types.hpp b/core/specfem/assembly/edge_types/dim2/edge_types.hpp index eee5f3ab2..ce9cb6a58 100644 --- a/core/specfem/assembly/edge_types/dim2/edge_types.hpp +++ b/core/specfem/assembly/edge_types/dim2/edge_types.hpp @@ -281,4 +281,14 @@ template <> struct edge_types { (EdgeViewType::HostMirror, h_coupled_edges))) }; +specfem::assembly::edge_types< + specfem::element::dimension_tag::dim2>::EdgeViewType::HostMirror +edge_view_from_collected_edges( + const std::string &label, + const std::vector< + specfem::mesh_entity::edge > + &collected_edges, + const specfem::mesh_entity::element + &element); + } // namespace specfem::assembly diff --git a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp index 9982c61de..036c3abf8 100644 --- a/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp +++ b/tests/unit-tests/algorithms/dim2/coupling_integral/dshape.cpp @@ -11,112 +11,6 @@ #include "utilities/include/fixture/nonconforming_interface/quadrature.hpp" #include -template -class ChunkEdgeIterator : public specfem::execution::TeamThreadRangePolicy< - Kokkos::TeamPolicy<>::member_type, int> { -private: - using KokkosIndexType = Kokkos::TeamPolicy<>::member_type; - -public: - using base_type = - specfem::execution::TeamThreadRangePolicy; - using execution_space = typename base_type::execution_space; - using index_type = - specfem::execution::EdgePointIndex; - - KOKKOS_INLINE_FUNCTION - ChunkEdgeIterator(const EdgeTypesView &edge_types, const int &nedges, - const int &ngllz, const int &ngllx, - const KokkosIndexType &team_member) - : edge_types(edge_types), nedges(nedges), ngllz(ngllz), ngllx(ngllz), - num_points(std::max(ngllx, ngllz)), - base_type(team_member, nedges * std::max(ngllx, ngllz)), - global_offset(0) {} - - KOKKOS_INLINE_FUNCTION - const index_type - operator()(const typename base_type::policy_index_type &i) const { - const int iedge = i % nedges; - const int ipoint = i / nedges; - - const specfem::mesh_entity::dim2::type edge_type = edge_types(iedge); - const bool is_leftright = - (edge_type == specfem::mesh_entity::dim2::type::left || - edge_type == specfem::mesh_entity::dim2::type::right); - const int num_points_norm = is_leftright ? ngllx : ngllz; - - const int inorm = (edge_type == specfem::mesh_entity::dim2::type::bottom || - edge_type == specfem::mesh_entity::dim2::type::left) - ? 0 - : num_points_norm - 1; - - const int iz = is_leftright ? ipoint : inorm; - const int ix = is_leftright ? inorm : ipoint; - - return index_type( - specfem::point::edge_index( - global_offset + iedge /*ispec*/, global_offset + iedge /*iedge*/, - ipoint, iz, ix, edge_type), - iedge, i); - } - - const int nedges; - -private: - EdgeTypesView edge_types; - const int global_offset; - const int ngllz; - const int ngllx; - const int num_points; -}; - -template class ChunkEdgeIndex { -public: - static constexpr auto accessor_type = - specfem::datatype::AccessorType::chunk_edge; - using KokkosIndexType = Kokkos::TeamPolicy<>::member_type; - using iterator_type = ChunkEdgeIterator; - - /** - * @brief Get Kokkos team member index. - * @return Reference to Kokkos team member - */ - KOKKOS_INLINE_FUNCTION - constexpr const KokkosIndexType &get_policy_index() const { - return this->kokkos_index; - } - - /** - * @brief Construct chunk edge index. - * @param nedges Number of edges in chunk - * @param kokkos_index Kokkos team member - */ - KOKKOS_INLINE_FUNCTION - ChunkEdgeIndex(const EdgeTypesView &edge_types, const int &nedges, - const int &ngllz, const int &ngllx, - const KokkosIndexType &kokkos_index) - : kokkos_index(kokkos_index), _nedges(nedges), - iterator(edge_types, nedges, ngllz, ngllx, kokkos_index) {} - - /** - * @brief Get number of edges. - * @return Edge count - */ - KOKKOS_INLINE_FUNCTION int nedges() const { return _nedges; } - -private: - int _nedges; ///< Number of edges in the chunk - KokkosIndexType kokkos_index; /**< Kokkos team member for this chunk */ - iterator_type iterator; - -public: - KOKKOS_INLINE_FUNCTION const iterator_type &get_iterator() const { - return iterator; - } -}; - /** * @brief Patches assembly::nonconforming_interfaces to not require mesh and * edge types. @@ -318,22 +212,16 @@ void execute_simple_dshape_test() { specfem::assembly::EdgeView edgelist( "dshape::edgelist", num_edges, nquad_element); { - specfem::assembly::EdgeView h_edgelist( - "dshape::h_edgelist", num_edges, nquad_element); + + std::vector > edge_collect; const auto element = specfem::mesh_entity::element(ngllz, ngllx); for (int iedge = 0; iedge < num_edges; iedge++) { - h_edgelist.element_index(iedge) = iedge; - h_edgelist.edge_index(iedge) = iedge; - h_edgelist.edge_types(iedge) = h_edge_types(iedge); - - for (int ipoint = 0; ipoint < nquad_element; ipoint++) { - const auto [iz, ix] = - element.map_coordinates(h_edge_types(iedge), ipoint); - h_edgelist.iz(iedge, ipoint) = iz; - h_edgelist.ix(iedge, ipoint) = ix; - } + edge_collect.push_back({ iedge, iedge, side, false }); } + const auto h_edgelist = specfem::assembly::edge_view_from_collected_edges( + "dshape::edgelist_host_mirror", edge_collect, element); + // specfem::assembly::edge_types::deep_copy(edgelist, // h_edgelist); Kokkos::deep_copy(edgelist.element_index, h_edgelist.element_index); @@ -445,10 +333,23 @@ void execute_simple_dshape_test() { edgelist, nonconforming_interfaces, ngllz, ngllx, lagrange_derivative.xi, intersection_contra_normal); - Kokkos::parallel_for( - "SimpleDShapeTest", Kokkos::TeamPolicy<>(num_edges, 1, 1), - KOKKOS_LAMBDA(const Kokkos::TeamPolicy<>::member_type &team_member) { - const int chunk_index = team_member.league_rank(); + using default_parallel_config = + specfem::parallel_configuration::default_chunk_edge_config< + dimension_tag, Kokkos::DefaultExecutionSpace>; + // override parallel config to have test chunk size. + using parallel_config = specfem::parallel_configuration::edge_chunk_config< + dimension_tag, chunk_size, default_parallel_config::execution_space>; + specfem::execution::ChunkedEdgeIterator chunk(parallel_config(), edgelist); + + specfem::execution::for_each_level( + "specfem::compute::shape_function_normal_derivative", chunk, + KOKKOS_LAMBDA( + const typename decltype(chunk)::index_type &chunk_iterator_index) { + const auto &iter_chunk_index = chunk_iterator_index.get_index(); + const auto &team = iter_chunk_index.get_policy_index(); + const int &num_edges = iter_chunk_index.nedges(); + + const int chunk_index = team.league_rank(); const FunctionType F( Kokkos::subview(function_view, Kokkos::make_pair(chunk_index * chunk_size, @@ -460,29 +361,10 @@ void execute_simple_dshape_test() { (chunk_index + 1) * chunk_size), Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); specfem::algorithms::coupling_integral_dnshape( - ngllz, ngllx, - ChunkEdgeIndex(Kokkos::subview(edge_types, - Kokkos::make_pair( - chunk_index * chunk_size, - (chunk_index + 1) * chunk_size)), - num_edges, ngllz, ngllx, team_member), - F, intersection_factor, SFND, + ngllz, ngllx, iter_chunk_index, F, intersection_factor, SFND, [&](const auto &index, const auto &point) { computed_integrals(index.ispec, index.iz, index.ix) = point(0); }); - // specfem::algorithms::coupling_integral_dnshape( - // nonconforming_interfaces, ngllz, ngllx, lagrange_derivative, - // ChunkEdgeIndex(Kokkos::subview(edge_types, - // Kokkos::make_pair( - // chunk_index * chunk_size, - // (chunk_index + 1) * - // chunk_size)), - // num_edges, ngllz, ngllx, team_member), - // F, intersection_factor, intersection_contra_normal, - // transfer_function_self_derivative, - // [&](const auto &index, const auto &point) { - // computed_integrals(index.ispec, index.iz, index.ix) = point(0); - // }); }); const auto h_computed_integrals = Kokkos::create_mirror_view_and_copy( From 72a7a387f1bdd0d0b1830123fd10473019390179 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 17 Feb 2026 13:23:09 -0500 Subject: [PATCH 14/15] missing pragma once --- .../fixture/nonconforming_interface/lagrange_derivative.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp index 6b29572be..6eecd5a1a 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp @@ -1,3 +1,5 @@ +#pragma once + #include "../impl/descriptions.hpp" #include "initializers.hpp" #include "specfem_setup.hpp" From 1be3df5ded39f00a2d02286634f2a579654f994c Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 17 Feb 2026 15:22:00 -0500 Subject: [PATCH 15/15] lagrange_derivative.hpp: specfem_setup -> specfem/setup --- .../fixture/nonconforming_interface/lagrange_derivative.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp index 6eecd5a1a..784963491 100644 --- a/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp +++ b/tests/unit-tests/utilities/include/fixture/nonconforming_interface/lagrange_derivative.hpp @@ -2,7 +2,7 @@ #include "../impl/descriptions.hpp" #include "initializers.hpp" -#include "specfem_setup.hpp" +#include "specfem/setup.hpp" namespace specfem::test_fixture {