From b4881418494dac0bb5e683efecabe96510ca40bd Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Tue, 2 Jun 2026 22:58:23 -0400 Subject: [PATCH 1/5] start 3d interface container, fwd declare ncis --- .../assembly/nonconforming_interfaces.hpp | 28 +-- .../dim3/impl/interface_container.hpp | 194 ++++++++++++++++++ .../dim3/impl/interface_container.tpp | 155 ++++++++++++++ .../dim3/nonconforming_interfaces.cpp | 20 ++ .../dim3/nonconforming_interfaces.hpp | 88 ++++++++ .../assembly/nonconforming_interfaces/fwd.hpp | 27 +++ 6 files changed, 485 insertions(+), 27 deletions(-) create mode 100644 core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp create mode 100644 core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp create mode 100644 core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp create mode 100644 core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp create mode 100644 core/specfem/assembly/nonconforming_interfaces/fwd.hpp diff --git a/core/specfem/assembly/nonconforming_interfaces.hpp b/core/specfem/assembly/nonconforming_interfaces.hpp index 8984312c30..aa4d4c5843 100644 --- a/core/specfem/assembly/nonconforming_interfaces.hpp +++ b/core/specfem/assembly/nonconforming_interfaces.hpp @@ -1,30 +1,4 @@ #pragma once -#include "specfem/element_connections.hpp" -#include "specfem/element_coupling.hpp" -#include "specfem/element_coupling/tags.hpp" -#include "specfem/enums.hpp" - -namespace specfem::assembly { - -namespace nonconforming_interfaces_impl { - -template -struct interface_container; - -} // namespace nonconforming_interfaces_impl - -/** - * @brief Information on coupled interfaces between two mediums - * @tparam DimensionTag Dimension of spectral elements - */ -template -struct nonconforming_interfaces; - -} // namespace specfem::assembly - #include "nonconforming_interfaces/dim2/nonconforming_interfaces.hpp" +#include "nonconforming_interfaces/dim3/nonconforming_interfaces.hpp" diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp new file mode 100644 index 0000000000..21bb90a8a8 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp @@ -0,0 +1,194 @@ +#pragma once + +#include "specfem/assembly/nonconforming_interfaces/fwd.hpp" + +#include "specfem/assembly/element_intersections.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/element/dimension.hpp" +#include "specfem/element_coupling/flux_scheme_configuration.hpp" +#include "specfem/element_coupling/tags.hpp" +#include "specfem/enums.hpp" +#include "specfem/execution.hpp" + +namespace specfem::assembly::nonconforming_interfaces_impl { + +template +struct interface_container; + +/** + * @brief Container for 3D nonconforming interface data storage and access + * + * Manages interface data between different physical media (elastic-acoustic) + * with specific boundary conditions. Stores edge factors and normal vectors + * for interface computations in 3D spectral element simulations. + * + * TODO: consider same physical media + * + * @tparam InterfaceTag Type of interface (ELASTIC_ACOUSTIC or ACOUSTIC_ELASTIC) + * @tparam BoundaryTag Boundary condition type (NONE, STACEY, etc.) + */ +template +struct interface_container< + specfem::element::dimension_tag::dim3, InterfaceTag, BoundaryTag, + specfem::element_connections::type::nonconforming, FluxSchemeTag> + : public specfem::data_access::Container< + specfem::data_access::ContainerType::face, + specfem::data_access::DataClassType::nonconforming_interface, + specfem::element::dimension_tag::dim3> { +public: + /** @brief Dimension tag for 2D specialization */ + constexpr static auto dimension_tag = specfem::element::dimension_tag::dim3; + /** @brief Interface type (elastic-acoustic or acoustic-elastic) */ + constexpr static auto interface_tag = InterfaceTag; + /** @brief Boundary condition type */ + constexpr static auto boundary_tag = BoundaryTag; + /** @brief Flux scheme type */ + constexpr static auto flux_scheme_tag = FluxSchemeTag; + /** @brief Medium type on the self side of the interface */ + constexpr static auto self_medium = + specfem::element_coupling::attributes::self_medium(); + /** @brief Medium type on the coupled side of the interface */ + constexpr static auto coupled_medium = + specfem::element_coupling::attributes::coupled_medium(); + +public: + /** @brief Base container type alias */ + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::face, + specfem::data_access::DataClassType::nonconforming_interface, + specfem::element::dimension_tag::dim3>; + /** @brief View type for edge scaling factors */ + using FaceFactorView = typename base_type::vector_type< + type_real, Kokkos::DefaultExecutionSpace::memory_space>; + /** @brief View type for edge normal vectors */ + using FaceNormalView = typename base_type::tensor_type< + type_real, Kokkos::DefaultExecutionSpace::memory_space>; + /** @brief View type for transfer function */ + using CoupledCoordinatesView = typename base_type::tensor_type< + type_real, Kokkos::DefaultExecutionSpace::memory_space>; + + /** @brief Device view for edge scaling factors */ + FaceFactorView face_factor; + /** @brief Device view for edge normal vectors */ + FaceNormalView face_normal; + /** @brief Device view for self nodes in coupled coordinates */ + CoupledCoordinatesView coupled_coordinates; + + /** @brief Host mirror for edge scaling factors */ + FaceFactorView::host_mirror_type h_face_factor; + /** @brief Host mirror for edge normal vectors */ + FaceNormalView::host_mirror_type h_face_normal; + /** @brief Device view for self nodes in coupled coordinates */ + CoupledCoordinatesView::host_mirror_type h_coupled_coordinates; + +public: + /** + * @brief Constructs interface container with mesh and geometry data + * + * @param ngllz Number of GLL points in z-direction + * @param ngllx Number of GLL points in x-direction + * @param element_intersections Element intersection information from mesh + * @param mesh Mesh connectivity and geometry + */ + interface_container( + const int &ngllz, const int &nglly, const int &ngllx, + const specfem::assembly::element_intersections< + specfem::element::dimension_tag::dim3> &element_intersections, + const specfem::assembly::mesh &mesh, + const specfem::element_coupling::flux_scheme_configuration + &flux_scheme_config = {}); + + /** @brief Default constructor */ + interface_container() = default; + + interface_container(const int &ngllz, const int &nglly, const int &ngllx, + const int &num_faces) + : face_factor("specfem::assembly::nonconforming_interfaces::face_factor", + num_faces, std::max(std::max(ngllz, nglly), ngllx), + std::max(std::max(ngllz, nglly), ngllx)), + h_face_factor(Kokkos::create_mirror_view(face_factor)), + face_normal("specfem::assembly::nonconforming_interfaces::face_normal", + num_faces, std::max(std::max(ngllz, nglly), ngllx), + std::max(std::max(ngllz, nglly), ngllx), + specfem::element::dimension::dim), + h_face_normal(Kokkos::create_mirror_view(face_normal)), + coupled_coordinates( + "specfem::assembly::nonconforming_interfaces::coupled_coordinates", + num_faces, std::max(std::max(ngllz, nglly), ngllx), + std::max(std::max(ngllz, nglly), ngllx), + specfem::element::dimension::dim - 1), + h_coupled_coordinates(Kokkos::create_mirror_view(coupled_coordinates)) { + }; + + /** + * @brief Loads interface data at specified index into point + * + * Template function that loads face factor and normal vector data + * from either device or host memory into the provided point object. + * + * @tparam on_device If true, loads from device memory; if false, from host + * @tparam IndexType Type of index (must have iface and ipoint members) + * @tparam PointType Type of point (must have face_factor and face_normal) + * @param index Face and point indices for data location + * @param point Output point object to store loaded data + */ + template + KOKKOS_FORCEINLINE_FUNCTION void + impl_load(const std::integral_constant< + specfem::datatype::AccessorType, + specfem::datatype::AccessorType::point> /* AccessorType */, + const IndexType &index, PointType &point) const { + + static_assert(specfem::data_access::is_point::value, + "impl_load only supports point accessors"); + + static_assert(specfem::data_access::is_point::value, + "impl_load requires point type for IndexType"); + + static_assert(specfem::data_access::is_face_index::value, + "impl_load requires face index type for IndexType"); + + static_assert( + specfem::data_access::is_conforming_interface::value, + "impl_load requires conforming interface point type for PointType"); + + if constexpr (on_device) { + point.face_factor = + face_factor(index.iface, index.ipoint_i, index.ipoint_j); + point.face_normal(0) = + face_normal(index.iface, index.ipoint_i, index.ipoint_j, 0); + point.face_normal(1) = + face_normal(index.iface, index.ipoint_i, index.ipoint_j, 1); + point.face_normal(2) = + face_normal(index.iface, index.ipoint_i, index.ipoint_j, 2); + point.coupled_coordinates(0) = + coupled_coordinates(index.iface, index.ipoint_i, index.ipoint_j, 0); + point.coupled_coordinates(1) = + coupled_coordinates(index.iface, index.ipoint_i, index.ipoint_j, 1); + } else { + point.face_factor = + h_face_factor(index.iface, index.ipoint_i, index.ipoint_j); + point.face_normal(0) = + h_face_normal(index.iface, index.ipoint_i, index.ipoint_j, 0); + point.face_normal(1) = + h_face_normal(index.iface, index.ipoint_i, index.ipoint_j, 1); + point.face_normal(2) = + h_face_normal(index.iface, index.ipoint_i, index.ipoint_j, 2); + point.coupled_coordinates(0) = + h_coupled_coordinates(index.iface, index.ipoint_i, index.ipoint_j, 0); + point.coupled_coordinates(1) = + h_coupled_coordinates(index.iface, index.ipoint_i, index.ipoint_j, 1); + } + return; + } +}; +} // namespace specfem::assembly::nonconforming_interfaces_impl diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp new file mode 100644 index 0000000000..0330ac9777 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp @@ -0,0 +1,155 @@ +#pragma once + +// #include "compute_intersection.hpp" +// #include "compute_intersection.tpp" +#include "interface_container.hpp" + +template +specfem::assembly::nonconforming_interfaces_impl::interface_container< + specfem::element::dimension_tag::dim3, InterfaceTag, BoundaryTag, + specfem::element_connections::type::nonconforming, FluxSchemeTag>:: + interface_container( + const int &ngllz, const int &nglly, const int &ngllx, + const specfem::assembly::element_intersections< + specfem::element::dimension_tag::dim3> &element_intersections, + const specfem::assembly::mesh &mesh, + const specfem::element_coupling::flux_scheme_configuration + &flux_scheme_config) { + + if (ngllz <= 0 || nglly <= 0 || ngllx <= 0) { + KOKKOS_ABORT_WITH_LOCATION("Invalid GLL grid size"); + } + + if (ngllz != ngllx || ngllz != nglly) { + KOKKOS_ABORT_WITH_LOCATION( + "The number of GLL points in z, y, and x must be the same."); + } + const int ngll = std::max(std::max(ngllz, nglly), ngllx); + constexpr int ndim = specfem::element::dimension::dim; + + const auto [self_faces, coupled_faces] = + element_intersections.get_intersections_on_host( + specfem::element_connections::type::nonconforming, InterfaceTag, + BoundaryTag, FluxSchemeTag); + + const auto &N = self_faces.N; + const auto &weights = mesh.h_weights; + + + const int num_faces = 0; // TODO count them + + face_factor = + FaceFactorView("specfem::assembly::nonconforming_interfaces::face_factor", + num_faces, ngll, ngll); + h_face_factor = Kokkos::create_mirror_view(face_factor); + + face_normal = + FaceNormalView("specfem::assembly::nonconforming_interfaces::face_normal", + num_faces, ngll, ngll, ndim); + h_face_normal = Kokkos::create_mirror_view(face_normal); + + coupled_coordinates = CoupledCoordinatesView( + "specfem::assembly::nonconforming_interfaces::coupled_coordinates", + num_faces, ngll, ngll, ndim - 1); + h_coupled_coordinates = Kokkos::create_mirror_view(coupled_coordinates); + + // // used when computing transfer functions + // const Kokkos::View *, + // Kokkos::HostSpace> + // icoorg("icoorg", mesh.ngnod); + // const Kokkos::View *, + // Kokkos::HostSpace> + // jcoorg("jcoorg", mesh.ngnod); + + // for (int i = 0; i < N; ++i) { + // const int ispec = self_edges(i).element_index; + // const auto iedge_type = self_edges(i).edge_type; + // const int jspec = coupled_edges(i).element_index; + // const auto jedge_type = coupled_edges(i).edge_type; + // for (int inod = 0; inod < mesh.ngnod; inod++) { + // icoorg(inod).x = mesh.h_control_node_coord(0, ispec, inod); + // icoorg(inod).z = mesh.h_control_node_coord(1, ispec, inod); + // jcoorg(inod).x = mesh.h_control_node_coord(0, jspec, inod); + // jcoorg(inod).z = mesh.h_control_node_coord(1, jspec, inod); + // } + // auto transfer_subview = + // Kokkos::subview(h_transfer_function, i, Kokkos::ALL, Kokkos::ALL); + // auto transfer_subview_other = + // Kokkos::subview(h_transfer_function_other, i, Kokkos::ALL, + // Kokkos::ALL); + // specfem::assembly::nonconforming_interfaces_impl::set_transfer_functions( + // icoorg, jcoorg, iedge_type, jedge_type, interface_knots, mesh.h_xi, + // transfer_subview, transfer_subview_other); + // // compute normal on edge + // const int npoints = element.number_of_points_on_orientation(iedge_type); + + // // compute factor by finding first derivative of position + // // along the edge and multiplying by the quadrature weight + // const Kokkos::View dr_intersection( + // "dr_intersection", nquad_intersection, 2); + + // for (int iquad = 0; iquad < nquad_intersection; iquad++) { + // dr_intersection(iquad, 0) = 0; + // dr_intersection(iquad, 1) = 0; + // } + // for (int iknot = 0; iknot < nquad_intersection; iknot++) { + // // get local coordinate (we can recover this from the transfer function + // by + // // interpolating x) + // type_real local_coord = 0; + // for (int ipoint = 0; ipoint < npoints; ipoint++) { + // local_coord += transfer_subview(iknot, ipoint) * mesh.h_xi(ipoint); + // } + + // // get global coordinate -- we interpolate against shape prime + // const auto [xi, gamma] = [&]() -> std::pair { + // if (iedge_type == specfem::mesh_entity::dim3::type::bottom) { + // return { local_coord, -1 }; + // } else if (iedge_type == specfem::mesh_entity::dim3::type::right) { + // return { 1, local_coord }; + // } else if (iedge_type == specfem::mesh_entity::dim3::type::top) { + // return { local_coord, 1 }; + // } else { + // return { -1, local_coord }; + // } + // }(); + // const auto loc = + // jacobian::compute_locations(icoorg, mesh.ngnod, xi, gamma); + + // // accumulate derivative at each quadrature point + // for (int iquad = 0; iquad < nquad_intersection; iquad++) { + // dr_intersection(iquad, 0) += interface_deriv(iquad, iknot) * loc.x; + // dr_intersection(iquad, 1) += interface_deriv(iquad, iknot) * loc.z; + // } + // } + + // // convert dr to ds and multiply by weights + // for (int iquad = 0; iquad < nquad_intersection; iquad++) { + // this->h_intersection_factor(i, iquad) = + // interface_weights(iquad) * + // std::sqrt(dr_intersection(iquad, 0) * dr_intersection(iquad, 0) + + // dr_intersection(iquad, 1) * dr_intersection(iquad, 1)); + + // type_real nx = dr_intersection(iquad, 1); + // type_real nz = -dr_intersection(iquad, 0); + // type_real mag = std::sqrt(nx * nx + nz * nz); + + // // previous computation was a 90 deg clockwise rotation. it should be + // CCW + // // for these cases: + // if (iedge_type == specfem::mesh_entity::dim3::type::top || + // iedge_type == specfem::mesh_entity::dim3::type::left) { + // mag *= -1; + // } + // this->h_intersection_normal(i, iquad, 0) = nx / mag; + // this->h_intersection_normal(i, iquad, 1) = nz / mag; + // } + // } + Kokkos::deep_copy(face_factor, h_face_factor); + Kokkos::deep_copy(face_normal, h_face_normal); + Kokkos::deep_copy(coupled_coordinates, h_coupled_coordinates); +} diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp new file mode 100644 index 0000000000..e07e731d91 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp @@ -0,0 +1,20 @@ + +#include "specfem/assembly/nonconforming_interfaces.hpp" +#include "impl/interface_container.tpp" +#include "specfem/assembly/element_intersections.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/enums.hpp" + +specfem::assembly::nonconforming_interfaces< + specfem::element::dimension_tag::dim3>:: + nonconforming_interfaces( + const int ngllz, const int ngllx, + const specfem::assembly::element_intersections< + specfem::element::dimension_tag::dim3> &element_intersections, + const specfem::assembly::mesh &mesh, + const specfem::element_coupling::flux_scheme_configuration + &flux_scheme_config) + : interface_container([&]() { + return InterfaceContainerTemplateType( + ngllz, ngllx, element_intersections, mesh, flux_scheme_config); + }){}; diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp new file mode 100644 index 0000000000..841245fd32 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp @@ -0,0 +1,88 @@ +#pragma once + +#include "impl/interface_container.hpp" +#include "specfem/assembly/element_intersections.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/data_access.hpp" +#include "specfem/element_coupling/flux_scheme_configuration.hpp" +#include "specfem/element_coupling/tags.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros/tag_dispatch.hpp" +#include "specfem/tag_dispatch.hpp" +#include +#include + +namespace specfem::assembly { + +/** + * @brief Information on coupled interfaces between two mediums + * @tparam DimensionTag Dimension of spectral elements + */ +template +struct nonconforming_interfaces; + +template <> +class nonconforming_interfaces + : public specfem::data_access::Container< + specfem::data_access::ContainerType::edge, + specfem::data_access::DataClassType::nonconforming_interface, + specfem::element::dimension_tag::dim3> { +public: + static constexpr auto dimension_tag = specfem::element::dimension_tag::dim3; + +protected: + template + using InterfaceContainerType = + specfem::assembly::nonconforming_interfaces_impl::interface_container< + dimension_tag, InterfaceTag, BoundaryTag, ConnectionTag, + FluxSchemeTag>; + + template + using InterfaceContainerTemplateType = + InterfaceContainerType; + + static constexpr auto combinations = + DIMENSION_SET(dim3) * CONNECTION_SET(nonconforming) * + INTERFACE_SET(elastic_acoustic, acoustic_elastic) * + BOUNDARY_SET(none, acoustic_free_surface, stacey, + composite_stacey_dirichlet) * + FLUX_SCHEME_SET(natural); + + specfem::tag_dispatch::TypedStorage + interface_container; + +public: + nonconforming_interfaces( + const int ngllz, const int ngllx, + const specfem::assembly::element_intersections + &element_intersections, + const specfem::assembly::mesh &mesh, + const specfem::element_coupling::flux_scheme_configuration + &flux_scheme_config = {}); + + nonconforming_interfaces() = default; + + template + KOKKOS_INLINE_FUNCTION const + InterfaceContainerType & + get_interface_container() const { + using TagsType = specfem::tags::Tags; + return interface_container.template get(); + } +}; + +} // namespace specfem::assembly + +#include "data_access/load.hpp" diff --git a/core/specfem/assembly/nonconforming_interfaces/fwd.hpp b/core/specfem/assembly/nonconforming_interfaces/fwd.hpp new file mode 100644 index 0000000000..5788abe802 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/fwd.hpp @@ -0,0 +1,27 @@ +#pragma once + +#include "specfem/element/tags.hpp" +#include "specfem/element_connections/tags.hpp" +#include "specfem/element_coupling/tags.hpp" + +namespace specfem::assembly { + +namespace nonconforming_interfaces_impl { + +template +struct interface_container; + +} // namespace nonconforming_interfaces_impl + +/** + * @brief Information on coupled interfaces between two mediums + * @tparam DimensionTag Dimension of spectral elements + */ +template +struct nonconforming_interfaces; + +} // namespace specfem::assembly From c94d8bb97513eb61c7192ed9de69efef510698f5 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Wed, 3 Jun 2026 19:41:42 -0400 Subject: [PATCH 2/5] container complete --- .../dim2/impl/interface_container.hpp | 7 +- .../dim3/data_access/load.hpp | 108 +++++++++ .../dim3/impl/interface_container.hpp | 2 + .../dim3/impl/interface_container.tpp | 206 +++++++++--------- .../dim3/nonconforming_interfaces.cpp | 3 +- .../dim3/nonconforming_interfaces.hpp | 2 +- 6 files changed, 221 insertions(+), 107 deletions(-) create mode 100644 core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.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 9167669955..3d08e52f64 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim2/impl/interface_container.hpp @@ -3,7 +3,7 @@ #include "specfem/assembly/element_intersections.hpp" #include "specfem/assembly/mesh.hpp" -#include "specfem/data_access.hpp" +#include "specfem/assembly/nonconforming_interfaces/fwd.hpp" #include "specfem/element_coupling/flux_scheme_configuration.hpp" #include "specfem/element_coupling/tags.hpp" #include "specfem/enums.hpp" @@ -169,8 +169,7 @@ struct interface_container< "impl_load requires edge_index type for IndexType"); constexpr int nquad_intersection = - std::tuple_element_t<0, - std::tuple >::n_quad_intersection; + std::tuple_element_t<0, std::tuple>::n_quad_intersection; const auto assign = [&](auto &point, auto &iquad, auto &index) { using PointType = std::decay_t; @@ -204,7 +203,7 @@ struct interface_container< "impl_load requires chunk_edge type for IndexType"); constexpr int nquad_intersection = - std::tuple_element_t<0, std::tuple >::n_quad_intersection; + std::tuple_element_t<0, std::tuple>::n_quad_intersection; const auto factor_subview = [&]() { if constexpr (on_device) { diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp new file mode 100644 index 0000000000..6bab9aefb6 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp @@ -0,0 +1,108 @@ +#pragma once + +#include "specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp" +#include "specfem/data_access.hpp" + +namespace specfem::assembly { + +/** + * @defgroup CoupledInterfaceDataAccess + * @brief Data access functions for coupled interface computation data + * + */ + +/** + * @brief Load interface data from container to point on host + * + * Loads coupled interface data using compile-time dispatch based on the point's + * template parameters (connection, interface, and boundary types). + * + * @ingroup CoupledInterfaceDataAccess + * + * @tparam IndexType Face index type + * @tparam ContainerType Coupled interfaces container type + * @tparam AccessorType Interface point (or face) type + * + * @param index Face index specifying the interface location + * @param container Coupled interfaces container holding interface data + * @param accessor Point or chunk_face object where loaded data will be stored + * + * @pre index refers to valid mesh face + * @pre container is properly initialized + * @pre point type matches supported interface combinations + * + * @note For host-side computations only. Use load_on_device for device code. + */ +template ::value) && + (specfem::data_access::is_nonconforming_interface< + ContainerType>::value)), + int> = 0> +inline void load_on_host(const IndexType &index, const ContainerType &container, + AccessorType &accessor) { + + static_assert( + specfem::data_access::CheckCompatibility::value, + "Incompatible types in load_on_host"); + + using accessor_dispatch = + std::integral_constant; + + container + .template get_interface_container() + .template impl_load(accessor_dispatch(), index, accessor); +} + +/** + * @brief Load interface data from container to face on device + * + * Loads coupled interface data using compile-time dispatch based on the face's + * template parameters (connection, interface, and boundary types). + * + * @ingroup CoupledInterfaceDataAccess + * + * @tparam IndexType Face index type + * @tparam ContainerType Coupled interfaces container type + * @tparam AccessorType Interface point (or face) type + * + * @param index Face index specifying the interface location + * @param container Coupled interfaces container holding interface data + * @param accessor Point or chunk_face object where loaded data will be stored + * + * @pre index refers to valid mesh face + * @pre container is properly initialized + * @pre point type matches supported interface combinations + * + * @note For device-side computations only. Use load_on_host for host code. + */ +template ::value) && + (specfem::data_access::is_nonconforming_interface< + ContainerType>::value)), + int> = 0> +KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index, + const ContainerType &container, + AccessorType &accessor) { + static_assert( + specfem::data_access::CheckCompatibility::value, + "Incompatible types in load_on_device"); + + using accessor_dispatch = + std::integral_constant; + container + .template get_interface_container() + .template impl_load(accessor_dispatch(), index, accessor); + + return; +} +} // namespace specfem::assembly diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp index 21bb90a8a8..ed711f535e 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp @@ -3,6 +3,7 @@ #include "specfem/assembly/nonconforming_interfaces/fwd.hpp" #include "specfem/assembly/element_intersections.hpp" +#include "specfem/assembly/jacobian_matrix.hpp" #include "specfem/assembly/mesh.hpp" #include "specfem/data_access/container.hpp" #include "specfem/element/dimension.hpp" @@ -103,6 +104,7 @@ struct interface_container< const int &ngllz, const int &nglly, const int &ngllx, const specfem::assembly::element_intersections< specfem::element::dimension_tag::dim3> &element_intersections, + const specfem::assembly::jacobian_matrix &jacobian_matrix, const specfem::assembly::mesh &mesh, const specfem::element_coupling::flux_scheme_configuration &flux_scheme_config = {}); diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp index 0330ac9777..ed6c9b2567 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp @@ -1,8 +1,9 @@ #pragma once -// #include "compute_intersection.hpp" -// #include "compute_intersection.tpp" #include "interface_container.hpp" +#include "specfem/algorithms/locate_point.hpp" +#include "specfem/point/global_coordinates.hpp" +#include template &element_intersections, + const specfem::assembly::jacobian_matrix + &jacobian_matrix, const specfem::assembly::mesh &mesh, const specfem::element_coupling::flux_scheme_configuration &flux_scheme_config) { @@ -34,12 +37,9 @@ specfem::assembly::nonconforming_interfaces_impl::interface_container< specfem::element_connections::type::nonconforming, InterfaceTag, BoundaryTag, FluxSchemeTag); - const auto &N = self_faces.N; + const auto &num_faces = self_faces.N; const auto &weights = mesh.h_weights; - - const int num_faces = 0; // TODO count them - face_factor = FaceFactorView("specfem::assembly::nonconforming_interfaces::face_factor", num_faces, ngll, ngll); @@ -55,100 +55,106 @@ specfem::assembly::nonconforming_interfaces_impl::interface_container< num_faces, ngll, ngll, ndim - 1); h_coupled_coordinates = Kokkos::create_mirror_view(coupled_coordinates); - // // used when computing transfer functions - // const Kokkos::View *, - // Kokkos::HostSpace> - // icoorg("icoorg", mesh.ngnod); - // const Kokkos::View *, - // Kokkos::HostSpace> - // jcoorg("jcoorg", mesh.ngnod); - - // for (int i = 0; i < N; ++i) { - // const int ispec = self_edges(i).element_index; - // const auto iedge_type = self_edges(i).edge_type; - // const int jspec = coupled_edges(i).element_index; - // const auto jedge_type = coupled_edges(i).edge_type; - // for (int inod = 0; inod < mesh.ngnod; inod++) { - // icoorg(inod).x = mesh.h_control_node_coord(0, ispec, inod); - // icoorg(inod).z = mesh.h_control_node_coord(1, ispec, inod); - // jcoorg(inod).x = mesh.h_control_node_coord(0, jspec, inod); - // jcoorg(inod).z = mesh.h_control_node_coord(1, jspec, inod); - // } - // auto transfer_subview = - // Kokkos::subview(h_transfer_function, i, Kokkos::ALL, Kokkos::ALL); - // auto transfer_subview_other = - // Kokkos::subview(h_transfer_function_other, i, Kokkos::ALL, - // Kokkos::ALL); - // specfem::assembly::nonconforming_interfaces_impl::set_transfer_functions( - // icoorg, jcoorg, iedge_type, jedge_type, interface_knots, mesh.h_xi, - // transfer_subview, transfer_subview_other); - // // compute normal on edge - // const int npoints = element.number_of_points_on_orientation(iedge_type); - - // // compute factor by finding first derivative of position - // // along the edge and multiplying by the quadrature weight - // const Kokkos::View dr_intersection( - // "dr_intersection", nquad_intersection, 2); - - // for (int iquad = 0; iquad < nquad_intersection; iquad++) { - // dr_intersection(iquad, 0) = 0; - // dr_intersection(iquad, 1) = 0; - // } - // for (int iknot = 0; iknot < nquad_intersection; iknot++) { - // // get local coordinate (we can recover this from the transfer function - // by - // // interpolating x) - // type_real local_coord = 0; - // for (int ipoint = 0; ipoint < npoints; ipoint++) { - // local_coord += transfer_subview(iknot, ipoint) * mesh.h_xi(ipoint); - // } - - // // get global coordinate -- we interpolate against shape prime - // const auto [xi, gamma] = [&]() -> std::pair { - // if (iedge_type == specfem::mesh_entity::dim3::type::bottom) { - // return { local_coord, -1 }; - // } else if (iedge_type == specfem::mesh_entity::dim3::type::right) { - // return { 1, local_coord }; - // } else if (iedge_type == specfem::mesh_entity::dim3::type::top) { - // return { local_coord, 1 }; - // } else { - // return { -1, local_coord }; - // } - // }(); - // const auto loc = - // jacobian::compute_locations(icoorg, mesh.ngnod, xi, gamma); - - // // accumulate derivative at each quadrature point - // for (int iquad = 0; iquad < nquad_intersection; iquad++) { - // dr_intersection(iquad, 0) += interface_deriv(iquad, iknot) * loc.x; - // dr_intersection(iquad, 1) += interface_deriv(iquad, iknot) * loc.z; - // } - // } - - // // convert dr to ds and multiply by weights - // for (int iquad = 0; iquad < nquad_intersection; iquad++) { - // this->h_intersection_factor(i, iquad) = - // interface_weights(iquad) * - // std::sqrt(dr_intersection(iquad, 0) * dr_intersection(iquad, 0) + - // dr_intersection(iquad, 1) * dr_intersection(iquad, 1)); - - // type_real nx = dr_intersection(iquad, 1); - // type_real nz = -dr_intersection(iquad, 0); - // type_real mag = std::sqrt(nx * nx + nz * nz); - - // // previous computation was a 90 deg clockwise rotation. it should be - // CCW - // // for these cases: - // if (iedge_type == specfem::mesh_entity::dim3::type::top || - // iedge_type == specfem::mesh_entity::dim3::type::left) { - // mag *= -1; - // } - // this->h_intersection_normal(i, iquad, 0) = nx / mag; - // this->h_intersection_normal(i, iquad, 1) = nz / mag; - // } - // } + for (int iface = 0; iface < num_faces; ++iface) { + const int ispec = self_faces(iface).element_index; + const auto iface_type = self_faces(iface).face_type; + const int jspec = coupled_faces(iface).element_index; + const auto jface_type = coupled_faces(iface).face_type; + + // face local indices + int ipoint_i, ipoint_j, ipoint_ortho; + + // max along local indices + int nglli, ngllj; + + const auto [ix, iy, iz] = [&]() -> std::tuple { + if (iface_type == specfem::mesh_entity::dim3::type::bottom) { + ipoint_ortho = 0; + nglli = ngllx; + ngllj = nglly; + return { ipoint_i, ipoint_j, ipoint_ortho }; + } else if (iface_type == specfem::mesh_entity::dim3::type::right) { + ipoint_ortho = ngllx - 1; + nglli = nglly; + ngllj = ngllz; + return { ipoint_ortho, ipoint_i, ipoint_j }; + } else if (iface_type == specfem::mesh_entity::dim3::type::top) { + ipoint_ortho = ngllz - 1; + nglli = ngllx; + ngllj = nglly; + return { ipoint_i, ipoint_j, ipoint_ortho }; + } else if (iface_type == specfem::mesh_entity::dim3::type::left) { + ipoint_ortho = 0; + nglli = nglly; + ngllj = ngllz; + return { ipoint_ortho, ipoint_i, ipoint_j }; + } else if (iface_type == specfem::mesh_entity::dim3::type::front) { + ipoint_ortho = 0; + nglli = ngllx; + ngllj = ngllz; + return { ipoint_i, ipoint_ortho, ipoint_j }; + } else { // back + ipoint_ortho = nglly - 1; + nglli = ngllx; + ngllj = ngllz; + return { ipoint_i, ipoint_ortho, ipoint_j }; + } + }(); + + for (int ipoint_i = 0; ipoint_i < nglli; ipoint_i++) { + for (int ipoint_j = 0; ipoint_j < ngllj; ipoint_j++) { + specfem::point::global_coordinates global_coord( + Kokkos::subview(mesh.h_coord, ispec, iz, iy, ix, Kokkos::ALL)); + + const auto [local_coords, found] = + specfem::algorithms::locate_point_impl::locate_point( + global_coord, mesh, jspec, jface_type); + + if (found) { + this->h_coupled_coordinates(iface, ipoint_i, ipoint_j, 0) = + local_coords.first; + this->h_coupled_coordinates(iface, ipoint_i, ipoint_j, 1) = + local_coords.second; + } else { + for (int idim = 0; idim < ndim - 1; idim++) { + this->h_coupled_coordinates(iface, ipoint_i, ipoint_j, idim) = NAN; + } + } + + specfem::point::jacobian_matrix + point_jacobian_matrix; + specfem::point::index + point_index{ ispec, iz, iy, ix }; + specfem::assembly::load_on_host(point_index, jacobian_matrix, + point_jacobian_matrix); + + const auto dn = point_jacobian_matrix.compute_normal(iface_type); + this->h_face_normal(iface, ipoint_i, ipoint_j, 0) = dn(0); + this->h_face_normal(iface, ipoint_i, ipoint_j, 1) = dn(1); + this->h_face_normal(iface, ipoint_i, ipoint_j, 2) = dn(2); + this->h_face_factor(iface, ipoint_i, ipoint_j) = [&]() { + switch (iface_type) { + case specfem::mesh_entity::dim3::type::left: + case specfem::mesh_entity::dim3::type::right: + // Face in (iy, iz) plane; integrate over iy and iz + return mesh.h_weights(iy) * mesh.h_weights(iz); + case specfem::mesh_entity::dim3::type::bottom: + case specfem::mesh_entity::dim3::type::top: + // Face in (ix, iy) plane; integrate over ix and iy + return mesh.h_weights(ix) * mesh.h_weights(iy); + case specfem::mesh_entity::dim3::type::front: + case specfem::mesh_entity::dim3::type::back: + // Face in (ix, iz) plane; integrate over ix and iz + return mesh.h_weights(ix) * mesh.h_weights(iz); + default: + KOKKOS_ABORT_WITH_LOCATION("Invalid face type"); + return static_cast(0.0); + } + }(); + } + } + } Kokkos::deep_copy(face_factor, h_face_factor); Kokkos::deep_copy(face_normal, h_face_normal); Kokkos::deep_copy(coupled_coordinates, h_coupled_coordinates); diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp index e07e731d91..5cdfa14fd3 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp @@ -1,6 +1,5 @@ -#include "specfem/assembly/nonconforming_interfaces.hpp" -#include "impl/interface_container.tpp" +#include "nonconforming_interfaces.hpp" #include "specfem/assembly/element_intersections.hpp" #include "specfem/assembly/mesh.hpp" #include "specfem/enums.hpp" diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp index 841245fd32..52b3bdec47 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp @@ -24,7 +24,7 @@ struct nonconforming_interfaces; template <> class nonconforming_interfaces : public specfem::data_access::Container< - specfem::data_access::ContainerType::edge, + specfem::data_access::ContainerType::face, specfem::data_access::DataClassType::nonconforming_interface, specfem::element::dimension_tag::dim3> { public: From 3aa32c7626bf275e6031180f25d7f2a284cda7a8 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Thu, 4 Jun 2026 16:34:15 -0400 Subject: [PATCH 3/5] temporary element intersection change for nci test --- .../dim3/element_intersections.cpp | 23 ++++++++--- .../specfem/compute/impl/compute_coupling.tpp | 39 +++++++++++-------- core/specfem/macros/interface_iterators.hpp | 8 +++- 3 files changed, 46 insertions(+), 24 deletions(-) diff --git a/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp b/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp index f9cf518e56..b4ac2f0365 100644 --- a/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp +++ b/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp @@ -1,5 +1,6 @@ #include "specfem/assembly/element_intersections.hpp" #include "specfem/assembly/element_types.hpp" +#include "specfem/element_connections/tags.hpp" #include "specfem/enums.hpp" #include "specfem/macros.hpp" #include "specfem/mesh.hpp" @@ -34,8 +35,8 @@ specfem::assembly::element_intersections< // Collect self/coupled face vectors for each tag combination struct CollectedFaces { - std::vector > self_collect; - std::vector > coupled_collect; + std::vector> self_collect; + std::vector> coupled_collect; }; specfem::tag_dispatch::Storage @@ -45,15 +46,25 @@ specfem::assembly::element_intersections< constexpr auto coupled_medium = specfem::element_coupling::attributes< TagsType::dimension_tag, TagsType::interface_tag>::coupled_medium(); - std::vector > self_collect; - std::vector > coupled_collect; + std::vector> self_collect; + std::vector> coupled_collect; const auto &graph = mesh.graph(); // Filter out corresponding connections auto filter = [&graph](const auto &edge) { + if constexpr (TagsType::connection_tag == + specfem::element_connections::type::weakly_conforming) { + return false; + } else if constexpr (TagsType::connection_tag == + specfem::element_connections::type:: + weakly_conforming) { + return true; + } return graph[edge].connection == TagsType::connection_tag; - }; + }; // TODO yell at Hanson if this is still in when the PR is reviewed. + // This is to test the nonconforming kernel on a conforming mesh + // prior to a NC mesher. // Create a filtered graph view const auto &nc_graph = boost::make_filtered_graph(graph, filter); @@ -164,7 +175,7 @@ specfem::assembly::element_intersections< specfem::assembly::face_view_from_collected_faces( const std::string &label, const std::vector< - specfem::mesh_entity::face > + specfem::mesh_entity::face> &collected_faces, const specfem::mesh_entity::element &element) { diff --git a/core/specfem/compute/impl/compute_coupling.tpp b/core/specfem/compute/impl/compute_coupling.tpp index 4572713b8f..4ef8b1a88b 100644 --- a/core/specfem/compute/impl/compute_coupling.tpp +++ b/core/specfem/compute/impl/compute_coupling.tpp @@ -5,6 +5,7 @@ #include "specfem/assembly/assembly.hpp" #include "specfem/boundary_conditions.hpp" #include "specfem/chunk_edge.hpp" +#include "specfem/element/tags.hpp" #include "specfem/element_connections.hpp" #include "specfem/enums.hpp" #include "specfem/execution.hpp" @@ -58,9 +59,12 @@ void compute_coupling_core_weakly_conforming( const auto num_points = assembly.mesh.element_grid.ngllx; - using parallel_config = + using parallel_config = std::conditional< + dimension_tag == specfem::element::dimension_tag::dim2, specfem::parallel_configuration::default_chunk_edge_config< - dimension_tag, Kokkos::DefaultExecutionSpace>; + dimension_tag, Kokkos::DefaultExecutionSpace>, + specfem::parallel_configuration::default_chunk_face_config< + dimension_tag, Kokkos::DefaultExecutionSpace>>; using CoupledFieldType = typename specfem::element_coupling::attributes< dimension_tag, interface_tag>::template coupled_field_t; @@ -139,9 +143,12 @@ void compute_coupling_core_nonconforming( const auto num_points = assembly.mesh.element_grid.ngllx; - using parallel_config = + using parallel_config = std::conditional< + dimension_tag == specfem::element::dimension_tag::dim2, specfem::parallel_configuration::default_chunk_edge_config< - dimension_tag, Kokkos::DefaultExecutionSpace>; + dimension_tag, Kokkos::DefaultExecutionSpace>, + specfem::parallel_configuration::default_chunk_face_config< + dimension_tag, Kokkos::DefaultExecutionSpace>>; // As written, field types cannot readily be defined in attributes. Define // them here. @@ -159,7 +166,7 @@ void compute_coupling_core_nonconforming( using_simd>, specfem::chunk_edge::acceleration >; + using_simd>>; using CouplingTermsPack = specfem::chunk_edge::coupling_terms_pack< dimension_tag, interface_tag, boundary_tag, flux_scheme_tag, @@ -172,7 +179,7 @@ void compute_coupling_core_nonconforming( type_real, dimension_tag, parallel_config::chunk_size, NQuad_intersection, specfem::element::attributes::components, using_simd, Kokkos::DefaultExecutionSpace::scratch_memory_space, - Kokkos::MemoryTraits >; + Kokkos::MemoryTraits>; specfem::execution::ChunkedIntersectionIterator chunk( parallel_config(), self_intersections, coupled_intersections); @@ -260,11 +267,11 @@ void compute_coupling( specfem::tag_dispatch::for_each( specfem::tag_dispatch::dimension_set{} * - CONNECTION_SET(weakly_conforming, nonconforming) * - INTERFACE_SET(elastic_acoustic, acoustic_elastic) * - BOUNDARY_SET(none, acoustic_free_surface, stacey, - composite_stacey_dirichlet) * - FLUX_SCHEME_SET(natural), + CONNECTION_SET(weakly_conforming, nonconforming) * + INTERFACE_SET(elastic_acoustic, acoustic_elastic) * + BOUNDARY_SET(none, acoustic_free_surface, stacey, + composite_stacey_dirichlet) * + FLUX_SCHEME_SET(natural), [&]() { constexpr auto self_medium = specfem::element_coupling::attributes< ElementTags::dimension_tag, @@ -272,12 +279,10 @@ void compute_coupling( if constexpr (self_medium == Tags::medium_tag) { compute_coupling_core< NGLL, NGLL, - specfem::tags::Tags >( + specfem::tags::Tags< + ElementTags::dimension_tag, ElementTags::connection_tag, + WavefieldType, ElementTags::interface_tag, + ElementTags::boundary_tag, ElementTags::flux_scheme_tag>>( assembly); } }); diff --git a/core/specfem/macros/interface_iterators.hpp b/core/specfem/macros/interface_iterators.hpp index 9c5233668a..f218de730e 100644 --- a/core/specfem/macros/interface_iterators.hpp +++ b/core/specfem/macros/interface_iterators.hpp @@ -208,9 +208,15 @@ INTERFACE_TAG_ACOUSTIC_ELASTIC, \ BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET, FLUX_SCHEME_TAG_NATURAL))( \ (DIMENSION_TAG_DIM3, CONNECTION_TAG_WEAKLY_CONFORMING, \ + INTERFACE_TAG_ELASTIC_ACOUSTIC, BOUNDARY_TAG_NONE, \ + FLUX_SCHEME_TAG_NATURAL))( \ + (DIMENSION_TAG_DIM3, CONNECTION_TAG_WEAKLY_CONFORMING, \ + INTERFACE_TAG_ACOUSTIC_ELASTIC, BOUNDARY_TAG_NONE, \ + FLUX_SCHEME_TAG_NATURAL))( \ + (DIMENSION_TAG_DIM3, CONNECTION_TAG_NONCONFORMING, \ INTERFACE_TAG_ELASTIC_ACOUSTIC, BOUNDARY_TAG_NONE, \ FLUX_SCHEME_TAG_NATURAL))((DIMENSION_TAG_DIM3, \ - CONNECTION_TAG_WEAKLY_CONFORMING, \ + CONNECTION_TAG_NONCONFORMING, \ INTERFACE_TAG_ACOUSTIC_ELASTIC, \ BOUNDARY_TAG_NONE, FLUX_SCHEME_TAG_NATURAL)) /// \endcond From ed3fd205c68f4b95b3c7711eb797283561075252 Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Fri, 5 Jun 2026 15:49:17 -0400 Subject: [PATCH 4/5] changes for 3d FS displacement test to run --- .../locate_point/dim3/locate_point_impl.cpp | 6 +- .../assembly/assembly/dim3/assembly.cpp | 10 +- .../assembly/assembly/dim3/assembly.hpp | 38 ++-- .../dim3/element_intersections.cpp | 2 +- .../impl/load_access_functions.hpp | 44 +++++ .../fields/data_access/load_on_device.hpp | 7 +- .../dim3/data_access/load.hpp | 12 +- .../dim3/impl/interface_container.hpp | 2 +- .../dim3/impl/interface_container.tpp | 8 +- .../dim3/nonconforming_interfaces.cpp | 8 +- .../dim3/nonconforming_interfaces.hpp | 3 +- .../specfem/compute/impl/compute_coupling.tpp | 172 ++++++++++++------ core/specfem/data_access/accessor.hpp | 19 ++ core/specfem/element_coupling/accessor.hpp | 27 ++- core/specfem/macros/interface_iterators.hpp | 8 +- .../point/nonconforming_interface3d.hpp | 29 ++- core/specfem/tag_dispatch/is_valid.hpp | 6 +- 17 files changed, 286 insertions(+), 115 deletions(-) diff --git a/core/specfem/algorithms/locate_point/dim3/locate_point_impl.cpp b/core/specfem/algorithms/locate_point/dim3/locate_point_impl.cpp index a7d66e85f1..8135e81f67 100644 --- a/core/specfem/algorithms/locate_point/dim3/locate_point_impl.cpp +++ b/core/specfem/algorithms/locate_point/dim3/locate_point_impl.cpp @@ -591,9 +591,9 @@ std::pair, bool> locate_point( Kokkos::HostSpace> coorg("coorg", mesh.ngnod); for (int i = 0; i < mesh.ngnod; i++) { - coorg(i).x = mesh.h_control_node_coordinates(0, ispec, i); - coorg(i).y = mesh.h_control_node_coordinates(1, ispec, i); - coorg(i).z = mesh.h_control_node_coordinates(2, ispec, i); + coorg(i).x = mesh.h_control_node_coordinates(ispec, i, 0); + coorg(i).y = mesh.h_control_node_coordinates(ispec, i, 1); + coorg(i).z = mesh.h_control_node_coordinates(ispec, i, 2); } // initial guess of 0 (center of edge) diff --git a/core/specfem/assembly/assembly/dim3/assembly.cpp b/core/specfem/assembly/assembly/dim3/assembly.cpp index 6d435aabe1..4193ddb9be 100644 --- a/core/specfem/assembly/assembly/dim3/assembly.cpp +++ b/core/specfem/assembly/assembly/dim3/assembly.cpp @@ -9,10 +9,10 @@ specfem::assembly::assembly::assembly( const specfem::mesh::mesh &mesh, const specfem::quadrature::quadratures &quadratures, - std::vector > > + std::vector>> &sources, const std::vector< - std::shared_ptr > > + std::shared_ptr>> &receivers, const std::vector &stypes, const type_real t0, const type_real dt, const int max_timesteps, const int max_sig_step, @@ -76,6 +76,12 @@ specfem::assembly::assembly::assembly( ngllz, nglly, ngllx, this->element_intersections, this->jacobian_matrix, this->mesh }; + + this->nonconforming_interfaces = { + ngllz, nglly, ngllx, this->element_intersections, this->jacobian_matrix, + this->mesh + }; + this->fields = { this->mesh, this->element_types, simulation }; // if (allocate_boundary_values) diff --git a/core/specfem/assembly/assembly/dim3/assembly.hpp b/core/specfem/assembly/assembly/dim3/assembly.hpp index cdb90eb4a2..3ec3cd316a 100644 --- a/core/specfem/assembly/assembly/dim3/assembly.hpp +++ b/core/specfem/assembly/assembly/dim3/assembly.hpp @@ -20,6 +20,8 @@ #include "specfem/receivers.hpp" #include "specfem/source.hpp" +#include "specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp" + namespace specfem::io { class reader; } @@ -125,6 +127,9 @@ template <> struct assembly { specfem::assembly::conforming_interfaces conforming_interfaces; ///< Conforming interfaces between different media ///< (e.g., fluid-solid interface) + specfem::assembly::nonconforming_interfaces + nonconforming_interfaces; ///< Nonconforming interfaces between different + ///< media (e.g., fluid-solid interface) specfem::assembly::boundary_values boundary_values; ///< Field @@ -166,23 +171,22 @@ template <> struct assembly { * assignment if exists) * @param flux_scheme_config Flux scheme rules for nonconforming interfaces */ - assembly( - const specfem::mesh::mesh &mesh, - const specfem::quadrature::quadratures &quadratures, - std::vector > > - &sources, - const std::vector< - std::shared_ptr > > - &receivers, - const std::vector &stypes, const type_real t0, - const type_real dt, const int max_timesteps, const int max_sig_step, - const int nsteps_between_samples, - const specfem::simulation::type simulation, - const bool allocate_boundary_values, - const std::shared_ptr &property_reader, - const specfem::element_coupling::flux_scheme_configuration - &flux_scheme_config = - specfem::element_coupling::flux_scheme_configuration()); + assembly(const specfem::mesh::mesh &mesh, + const specfem::quadrature::quadratures &quadratures, + std::vector>> + &sources, + const std::vector< + std::shared_ptr>> + &receivers, + const std::vector &stypes, + const type_real t0, const type_real dt, const int max_timesteps, + const int max_sig_step, const int nsteps_between_samples, + const specfem::simulation::type simulation, + const bool allocate_boundary_values, + const std::shared_ptr &property_reader, + const specfem::element_coupling::flux_scheme_configuration + &flux_scheme_config = + specfem::element_coupling::flux_scheme_configuration()); assembly() = default; diff --git a/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp b/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp index b4ac2f0365..e7b0a99ed3 100644 --- a/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp +++ b/core/specfem/assembly/element_intersections/dim3/element_intersections.cpp @@ -58,7 +58,7 @@ specfem::assembly::element_intersections< return false; } else if constexpr (TagsType::connection_tag == specfem::element_connections::type:: - weakly_conforming) { + nonconforming) { return true; } return graph[edge].connection == TagsType::connection_tag; diff --git a/core/specfem/assembly/fields/data_access/impl/load_access_functions.hpp b/core/specfem/assembly/fields/data_access/impl/load_access_functions.hpp index 1d1df3afdd..552a656098 100644 --- a/core/specfem/assembly/fields/data_access/impl/load_access_functions.hpp +++ b/core/specfem/assembly/fields/data_access/impl/load_access_functions.hpp @@ -256,4 +256,48 @@ load_after_simd_dispatch(const std::false_type, const IndexType &index, return; } +template ::value && + specfem::data_access::is_chunk_face::value && + (specfem::data_access::is_field::value && ...)), + int> = 0> +KOKKOS_FORCEINLINE_FUNCTION void +load_after_simd_dispatch(const std::false_type, const IndexType &index, + const ContainerType &field, + AccessorTypes &...accessors) { + static_assert(!IndexType::using_simd, + "IndexType must not use SIMD in this overload"); + constexpr static auto MediumTag = + std::tuple_element_t<0, std::tuple>::medium_tag; + // Check that all accessors have the same medium tag + specfem::assembly::fields_impl::check_accessor_compatibility< + AccessorTypes...>(); + const auto ¤t_field = field.template get_field(); + + constexpr static int ncomponents = specfem::element::attributes< + std::tuple_element_t<0, std::tuple>::dimension_tag, + std::tuple_element_t<0, std::tuple>::medium_tag>:: + components; + + specfem::execution::for_each_level( + index.get_iterator(), + [&](const typename IndexType::iterator_type::index_type &iterator_index) { + const auto point_index = iterator_index.get_local_index(); + const int iglob = + field.template get_iglob(point_index); + for (int icomp = 0; icomp < ncomponents; ++icomp) { + (specfem::assembly::fields_impl::base_load_accessor< + on_device, AccessorTypes::data_class>( + iglob, icomp, current_field, + accessors(point_index.iface, point_index.ipoint_i, + point_index.ipoint_j, icomp)), + ...); + } + }); + + return; +} + } // namespace specfem::assembly::simulation_field_impl diff --git a/core/specfem/assembly/fields/data_access/load_on_device.hpp b/core/specfem/assembly/fields/data_access/load_on_device.hpp index 7557cdabff..09ccd886a5 100644 --- a/core/specfem/assembly/fields/data_access/load_on_device.hpp +++ b/core/specfem/assembly/fields/data_access/load_on_device.hpp @@ -22,7 +22,7 @@ KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index, AccessorTypes &...accessors) { constexpr static auto MediumTag = - std::tuple_element_t<0, std::tuple >::medium_tag; + std::tuple_element_t<0, std::tuple>::medium_tag; // Check that all accessors have the same medium tag fields_impl::check_accessor_compatibility(); @@ -40,7 +40,8 @@ template ::value) && (specfem::data_access::is_point::value || specfem::data_access::is_chunk_element::value || - specfem::data_access::is_chunk_edge::value) && + specfem::data_access::is_chunk_edge::value || + specfem::data_access::is_chunk_face::value) && (specfem::data_access::is_field::value && ...)), int> = 0> KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index, @@ -48,7 +49,7 @@ KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index, AccessorTypes &...accessors) { constexpr static auto MediumTag = - std::tuple_element_t<0, std::tuple >::medium_tag; + std::tuple_element_t<0, std::tuple>::medium_tag; // Check that all accessors have the same medium tag fields_impl::check_accessor_compatibility(); diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp index 6bab9aefb6..c25e708728 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/data_access/load.hpp @@ -52,9 +52,9 @@ inline void load_on_host(const IndexType &index, const ContainerType &container, IndexType::accessor_type>; container - .template get_interface_container() + .template get_interface_container< + AccessorType::interface_tag, AccessorType::boundary_tag, + AccessorType::connection_tag, AccessorType::flux_scheme_tag>() .template impl_load(accessor_dispatch(), index, accessor); } @@ -98,9 +98,9 @@ KOKKOS_FORCEINLINE_FUNCTION void load_on_device(const IndexType &index, std::integral_constant; container - .template get_interface_container() + .template get_interface_container< + AccessorType::interface_tag, AccessorType::boundary_tag, + AccessorType::connection_tag, AccessorType::flux_scheme_tag>() .template impl_load(accessor_dispatch(), index, accessor); return; diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp index ed711f535e..bf8fc5519e 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp @@ -160,7 +160,7 @@ struct interface_container< "impl_load requires face index type for IndexType"); static_assert( - specfem::data_access::is_conforming_interface::value, + specfem::data_access::is_nonconforming_interface::value, "impl_load requires conforming interface point type for PointType"); if constexpr (on_device) { diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp index ed6c9b2567..a077b2ede7 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp @@ -101,10 +101,12 @@ specfem::assembly::nonconforming_interfaces_impl::interface_container< } }(); - for (int ipoint_i = 0; ipoint_i < nglli; ipoint_i++) { - for (int ipoint_j = 0; ipoint_j < ngllj; ipoint_j++) { + for (ipoint_i = 0; ipoint_i < nglli; ipoint_i++) { + for (ipoint_j = 0; ipoint_j < ngllj; ipoint_j++) { specfem::point::global_coordinates global_coord( - Kokkos::subview(mesh.h_coord, ispec, iz, iy, ix, Kokkos::ALL)); + mesh.h_coord(ispec, iz, iy, ix, 0), + mesh.h_coord(ispec, iz, iy, ix, 1), + mesh.h_coord(ispec, iz, iy, ix, 2)); const auto [local_coords, found] = specfem::algorithms::locate_point_impl::locate_point( diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp index 5cdfa14fd3..d0d41372a7 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp @@ -1,5 +1,6 @@ #include "nonconforming_interfaces.hpp" +#include "impl/interface_container.tpp" #include "specfem/assembly/element_intersections.hpp" #include "specfem/assembly/mesh.hpp" #include "specfem/enums.hpp" @@ -7,13 +8,16 @@ specfem::assembly::nonconforming_interfaces< specfem::element::dimension_tag::dim3>:: nonconforming_interfaces( - const int ngllz, const int ngllx, + const int &ngllz, const int &nglly, const int &ngllx, const specfem::assembly::element_intersections< specfem::element::dimension_tag::dim3> &element_intersections, + const specfem::assembly::jacobian_matrix + &jacobian_matrix, const specfem::assembly::mesh &mesh, const specfem::element_coupling::flux_scheme_configuration &flux_scheme_config) : interface_container([&]() { return InterfaceContainerTemplateType( - ngllz, ngllx, element_intersections, mesh, flux_scheme_config); + ngllz, nglly, ngllx, element_intersections, jacobian_matrix, mesh, + flux_scheme_config); }){}; diff --git a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp index 52b3bdec47..23328e2f39 100644 --- a/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp @@ -59,9 +59,10 @@ class nonconforming_interfaces public: nonconforming_interfaces( - const int ngllz, const int ngllx, + const int &ngllz, const int &nglly, const int &ngllx, const specfem::assembly::element_intersections &element_intersections, + const specfem::assembly::jacobian_matrix &jacobian_matrix, const specfem::assembly::mesh &mesh, const specfem::element_coupling::flux_scheme_configuration &flux_scheme_config = {}); diff --git a/core/specfem/compute/impl/compute_coupling.tpp b/core/specfem/compute/impl/compute_coupling.tpp index ecca8b42c9..336c0bf22a 100644 --- a/core/specfem/compute/impl/compute_coupling.tpp +++ b/core/specfem/compute/impl/compute_coupling.tpp @@ -2,13 +2,18 @@ #include "compute_coupling.hpp" #include "specfem/algorithms.hpp" +#include "specfem/algorithms/transfer_interpolate.hpp" #include "specfem/assembly/assembly.hpp" #include "specfem/boundary_conditions.hpp" #include "specfem/chunk_edge.hpp" +#include "specfem/chunk_face.hpp" +#include "specfem/data_access/accessor.hpp" +#include "specfem/data_access/accessor/point_accessor.hpp" #include "specfem/element/tags.hpp" #include "specfem/element_connections.hpp" #include "specfem/enums.hpp" #include "specfem/execution.hpp" +#include "specfem/execution/for_each_level.hpp" #include "specfem/medium_physics.hpp" #include "specfem/parallel_configuration.hpp" #include "specfem/point.hpp" @@ -61,7 +66,7 @@ void compute_coupling_core_weakly_conforming( const auto num_points = assembly.mesh.element_grid.ngllx; - using parallel_config = std::conditional< + using parallel_config = std::conditional_t< dimension_tag == specfem::element::dimension_tag::dim2, specfem::parallel_configuration::default_chunk_edge_config< dimension_tag, Kokkos::DefaultExecutionSpace>, @@ -134,6 +139,8 @@ void compute_coupling_core_nonconforming( "Currently, we are enforcing only one flux scheme: natural"); const auto &nonconforming_interfaces = assembly.nonconforming_interfaces; + const auto &boundaries = assembly.boundaries; + const auto [self_intersections, coupled_intersections] = assembly.element_intersections.get_intersections_on_device( connection_tag, interface_tag, boundary_tag, flux_scheme_tag); @@ -145,7 +152,7 @@ void compute_coupling_core_nonconforming( const auto num_points = assembly.mesh.element_grid.ngllx; - using parallel_config = std::conditional< + using parallel_config = std::conditional_t< dimension_tag == specfem::element::dimension_tag::dim2, specfem::parallel_configuration::default_chunk_edge_config< dimension_tag, Kokkos::DefaultExecutionSpace>, @@ -161,29 +168,52 @@ void compute_coupling_core_nonconforming( specfem::element_coupling::attributes::coupled_medium(); using CoupledFieldType = std::conditional_t< - interface_tag == - specfem::element_coupling::interface_tag::acoustic_elastic, - specfem::chunk_edge::displacement, - specfem::chunk_edge::acceleration>; + dimension_tag == specfem::element::dimension_tag::dim2, + std::conditional_t< + interface_tag == + specfem::element_coupling::interface_tag::acoustic_elastic, + specfem::chunk_edge::displacement, + specfem::chunk_edge::acceleration>, + std::conditional_t< + interface_tag == + specfem::element_coupling::interface_tag::acoustic_elastic, + specfem::chunk_face::displacement, + specfem::chunk_face::acceleration>>; using CouplingTermsPack = specfem::element_coupling::accessor::coupling_terms_pack< dimension_tag, interface_tag, boundary_tag, flux_scheme_tag, parallel_config::chunk_size, NGLL, NQuad_intersection>; - using IntegrationFactor = + + // should the nonconforming transfer be computed self-pointwise, or is it an + // intersection-type, where the entire intersection needs to be computed + // together? + constexpr bool is_pointwise_coupling = + specfem::data_access::is_point::value; + + // teamwise integration data, only needed if not pointwise + using IntegrationFactor = std::conditional_t< + is_pointwise_coupling, specfem::data_access::EmptyAccessor, specfem::element_coupling::accessor::intersection_factor< dimension_tag, interface_tag, boundary_tag, flux_scheme_tag, - parallel_config::chunk_size, NQuad_intersection>; + parallel_config::chunk_size, NQuad_intersection>>; - using InterfaceFieldViewType = specfem::datatype::VectorChunkEdgeViewType< - type_real, dimension_tag, parallel_config::chunk_size, NQuad_intersection, - specfem::element::attributes::components, - using_simd, Kokkos::DefaultExecutionSpace::scratch_memory_space, - Kokkos::MemoryTraits>; + using InterfaceFieldViewType = std::conditional_t< + is_pointwise_coupling, specfem::data_access::EmptyAccessor, + specfem::datatype::VectorChunkEdgeViewType< + type_real, dimension_tag, parallel_config::chunk_size, + NQuad_intersection, + specfem::element::attributes::components, + using_simd, Kokkos::DefaultExecutionSpace::scratch_memory_space, + Kokkos::MemoryTraits>>; specfem::execution::ChunkedIntersectionIterator chunk( parallel_config(), self_intersections, coupled_intersections); @@ -210,40 +240,80 @@ void compute_coupling_core_nonconforming( specfem::assembly::load_on_device(coupled_chunk_index, field, coupled_field); - CouplingTermsPack interface_data(team); - - specfem::assembly::load_on_device( - self_chunk_index, nonconforming_interfaces, interface_data); - InterfaceFieldViewType interface_field(team.team_scratch(0)); - - team.team_barrier(); - specfem::medium_physics::compute_coupling( - self_chunk_index, interface_data, coupled_field, interface_field); - - IntegrationFactor integration_factor(team); - - specfem::assembly::load_on_device( - self_chunk_index, nonconforming_interfaces, integration_factor); - - team.team_barrier(); - - specfem::algorithms::coupling_integral( - assembly.nonconforming_interfaces, self_chunk_index, - interface_field, integration_factor, - [&](const auto &self_index, auto &self_field) { - specfem::point::boundary - point_boundary; - specfem::assembly::load_on_device(self_index, assembly.boundaries, - point_boundary); - if constexpr (boundary_tag == specfem::element::boundary_tag:: - acoustic_free_surface) { - specfem::boundary_conditions::apply_boundary_conditions( - point_boundary, self_field); - } - - specfem::assembly::atomic_add_on_device(self_index, field, - self_field); - }); + if constexpr (is_pointwise_coupling) { + specfem::execution::for_each_level( + chunk_index.get_iterator(), + [&](const typename std::decay_t:: + iterator_type::index_type &iterator_index) { + const auto &index = iterator_index.get_index(); + + using SelfFieldType = + specfem::point::acceleration>; + + CouplingTermsPack point_interface_data; + specfem::assembly::load_on_device(index.self_index, + nonconforming_interfaces, + point_interface_data); + // TEMPORARY until we get rid of non-static interpolants + point_interface_data.set_interpolants( + specfem::algorithms::LagrangeInterpolant(assembly.mesh.xi)); + + SelfFieldType self_field; + specfem::medium_physics::compute_coupling( + index.self_index, point_interface_data, coupled_field, + self_field); + + specfem::point::boundary + point_boundary; + specfem::assembly::load_on_device(index.self_index, boundaries, + point_boundary); + if constexpr (boundary_tag == specfem::element::boundary_tag:: + acoustic_free_surface) { + specfem::boundary_conditions::apply_boundary_conditions( + point_boundary, self_field); + } + specfem::assembly::atomic_add_on_device(index.self_index, field, + self_field); + }); + + } else { + + CouplingTermsPack interface_data(team); + + specfem::assembly::load_on_device( + self_chunk_index, nonconforming_interfaces, interface_data); + InterfaceFieldViewType interface_field(team.team_scratch(0)); + + team.team_barrier(); + specfem::medium_physics::compute_coupling( + self_chunk_index, interface_data, coupled_field, interface_field); + + IntegrationFactor integration_factor(team); + + specfem::assembly::load_on_device( + self_chunk_index, nonconforming_interfaces, integration_factor); + + team.team_barrier(); + + specfem::algorithms::coupling_integral( + assembly.nonconforming_interfaces, self_chunk_index, + interface_field, integration_factor, + [&](const auto &self_index, auto &self_field) { + specfem::point::boundary + point_boundary; + specfem::assembly::load_on_device( + self_index, assembly.boundaries, point_boundary); + if constexpr (boundary_tag == specfem::element::boundary_tag:: + acoustic_free_surface) { + specfem::boundary_conditions::apply_boundary_conditions( + point_boundary, self_field); + } + + specfem::assembly::atomic_add_on_device(self_index, field, + self_field); + }); + } }); return; diff --git a/core/specfem/data_access/accessor.hpp b/core/specfem/data_access/accessor.hpp index 99bb621646..e07c0954d2 100644 --- a/core/specfem/data_access/accessor.hpp +++ b/core/specfem/data_access/accessor.hpp @@ -54,6 +54,25 @@ struct is_codim1_chunk< T::dimension_tag == specfem::element::dimension_tag::dim3)>> : std::true_type {}; +/** + * @brief A data accessor that holds no data, and contains null shmem info, etc. + * + * This container can be used as a placeholder when one kernel needs some data + * type (say in a scratch view) only for certain template arguments. + */ +struct EmptyAccessor { + + /** + * @brief Capture any constructer configuration. Should do nothing. + */ + template EmptyAccessor(Args...){}; + /** + * @brief Get shared memory size requirement + * @return Size in bytes needed for scratch memory + */ + constexpr static int shmem_size() { return 0; } +}; + } // namespace specfem::data_access #include "accessor/chunk_edge.hpp" diff --git a/core/specfem/element_coupling/accessor.hpp b/core/specfem/element_coupling/accessor.hpp index 9340d10692..9da2f087c6 100644 --- a/core/specfem/element_coupling/accessor.hpp +++ b/core/specfem/element_coupling/accessor.hpp @@ -1,7 +1,8 @@ #pragma once #include "specfem/chunk_edge/nonconforming_interface.hpp" -#include "specfem/chunk_face/nonconforming_interface.hpp" +#include "specfem/data_access/accessor.hpp" +#include "specfem/point/nonconforming_interface3d.hpp" #include #include @@ -61,13 +62,8 @@ struct coupling_terms_pack< DimensionTag, InterfaceTag, BoundaryTag, FluxSchemeTag, NumberElements, NQuadIntersection, NQuadElement, std::tuple, std::enable_if_t> { - using type = specfem::chunk_face::NonconformingAccessorPack< - specfem::chunk_face::coupled_coordinates< - DimensionTag, NumberElements, NQuadElement, InterfaceTag, BoundaryTag, - FluxSchemeTag, KokkosViewArgs...>, - specfem::chunk_face::intersection_normal< - DimensionTag, InterfaceTag, BoundaryTag, FluxSchemeTag, - NumberElements, NQuadElement, KokkosViewArgs...>>; + using type = specfem::point::nonconforming_interface< + DimensionTag, NQuadElement, InterfaceTag, BoundaryTag, FluxSchemeTag>; }; // =========================================================================== @@ -91,6 +87,21 @@ struct intersection_factor< NQuadIntersection>; }; +template +struct intersection_factor< + DimensionTag, InterfaceTag, BoundaryTag, FluxSchemeTag, NumberElements, + NQuadIntersection, std::tuple, + std::enable_if_t> { + static_assert(sizeof...(KokkosViewArgs) == 0, + "This coupling_terms_pack does not have Kokkos-view-argument " + "passthrough implemented!"); + using type = specfem::data_access::EmptyAccessor; +}; + // =========================================================================== } // namespace impl diff --git a/core/specfem/macros/interface_iterators.hpp b/core/specfem/macros/interface_iterators.hpp index f218de730e..9c5233668a 100644 --- a/core/specfem/macros/interface_iterators.hpp +++ b/core/specfem/macros/interface_iterators.hpp @@ -208,15 +208,9 @@ INTERFACE_TAG_ACOUSTIC_ELASTIC, \ BOUNDARY_TAG_COMPOSITE_STACEY_DIRICHLET, FLUX_SCHEME_TAG_NATURAL))( \ (DIMENSION_TAG_DIM3, CONNECTION_TAG_WEAKLY_CONFORMING, \ - INTERFACE_TAG_ELASTIC_ACOUSTIC, BOUNDARY_TAG_NONE, \ - FLUX_SCHEME_TAG_NATURAL))( \ - (DIMENSION_TAG_DIM3, CONNECTION_TAG_WEAKLY_CONFORMING, \ - INTERFACE_TAG_ACOUSTIC_ELASTIC, BOUNDARY_TAG_NONE, \ - FLUX_SCHEME_TAG_NATURAL))( \ - (DIMENSION_TAG_DIM3, CONNECTION_TAG_NONCONFORMING, \ INTERFACE_TAG_ELASTIC_ACOUSTIC, BOUNDARY_TAG_NONE, \ FLUX_SCHEME_TAG_NATURAL))((DIMENSION_TAG_DIM3, \ - CONNECTION_TAG_NONCONFORMING, \ + CONNECTION_TAG_WEAKLY_CONFORMING, \ INTERFACE_TAG_ACOUSTIC_ELASTIC, \ BOUNDARY_TAG_NONE, FLUX_SCHEME_TAG_NATURAL)) /// \endcond diff --git a/core/specfem/point/nonconforming_interface3d.hpp b/core/specfem/point/nonconforming_interface3d.hpp index aed93f023b..7c011ef14c 100644 --- a/core/specfem/point/nonconforming_interface3d.hpp +++ b/core/specfem/point/nonconforming_interface3d.hpp @@ -56,14 +56,13 @@ struct nonconforming_interface interpolants; + + /** + * @brief Populates the interpolants. + */ template - KOKKOS_INLINE_FUNCTION nonconforming_interface( - const vector_type &coupled_coordinates, - const scalar_type &face_factor, - const vector_type &face_normal, - const LagrangeInterpolantType &lagrange_interpolant) - : coupled_coordinates(coupled_coordinates), face_factor(face_factor), - face_normal(face_normal) { + KOKKOS_INLINE_FUNCTION void + set_interpolants(const LagrangeInterpolantType &lagrange_interpolant) { // populate interpolants array. for (int igll = 0; igll < NGLL; igll++) { for (int idim = 0; idim < ndim - 1; idim++) { @@ -72,6 +71,16 @@ struct nonconforming_interface + KOKKOS_INLINE_FUNCTION nonconforming_interface( + const vector_type &coupled_coordinates, + const scalar_type &face_factor, + const vector_type &face_normal, + const LagrangeInterpolantType &lagrange_interpolant) + : coupled_coordinates(coupled_coordinates), face_factor(face_factor), + face_normal(face_normal) { + set_interpolants(lagrange_interpolant); + } // ========================= END TEMPORARY ========================= @@ -87,6 +96,12 @@ struct nonconforming_interface, TagTypes...> * @return The stored enum value. */ template constexpr auto get() const { - using T = typename std::tuple_element >::type; + using T = typename std::tuple_element>::type; return static_cast &>(*this).v; } @@ -341,9 +341,9 @@ constexpr bool is_valid_edge(EdgeTuple t) { if (i != I::elastic_acoustic && i != I::acoustic_elastic) return false; - // dim3: only weakly_conforming + none + // dim3: only none if (d == D::dim3) - return c == C::weakly_conforming && b == B::none; + return b == B::none; // dim2: boundary rules depend on interface direction if (i == I::elastic_acoustic) From 2c32a476157888689ebf40f0b06f71326b85889f Mon Sep 17 00:00:00 2001 From: Kentaro Hanson Date: Fri, 5 Jun 2026 18:12:46 -0400 Subject: [PATCH 5/5] rewrote outer parfor due to --- .../specfem/compute/impl/compute_coupling.tpp | 31 +++++++++++++------ 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/core/specfem/compute/impl/compute_coupling.tpp b/core/specfem/compute/impl/compute_coupling.tpp index 336c0bf22a..60a78b6e3f 100644 --- a/core/specfem/compute/impl/compute_coupling.tpp +++ b/core/specfem/compute/impl/compute_coupling.tpp @@ -240,40 +240,51 @@ void compute_coupling_core_nonconforming( specfem::assembly::load_on_device(coupled_chunk_index, field, coupled_field); - if constexpr (is_pointwise_coupling) { + + + + const auto &nonconforming_interfaces = + assembly.nonconforming_interfaces; + const auto& assembly_mesh_xi = assembly.mesh.xi; + const auto& boundaries = assembly.boundaries; + + + // internal to lambda, since nvcc could not compile otherwise + constexpr bool is_pointwise_coupling_ = is_pointwise_coupling; + if constexpr (is_pointwise_coupling_) { specfem::execution::for_each_level( - chunk_index.get_iterator(), - [&](const typename std::decay_t:: + self_chunk_index.get_iterator(), + [&](const typename std::decay_t:: iterator_type::index_type &iterator_index) { const auto &index = iterator_index.get_index(); + const auto &local_index = iterator_index.get_local_index(); using SelfFieldType = specfem::point::acceleration>; CouplingTermsPack point_interface_data; - specfem::assembly::load_on_device(index.self_index, - nonconforming_interfaces, - point_interface_data); + specfem::assembly::load_on_device( + index, nonconforming_interfaces, point_interface_data); // TEMPORARY until we get rid of non-static interpolants point_interface_data.set_interpolants( - specfem::algorithms::LagrangeInterpolant(assembly.mesh.xi)); + specfem::algorithms::LagrangeInterpolant(assembly_mesh_xi)); SelfFieldType self_field; specfem::medium_physics::compute_coupling( - index.self_index, point_interface_data, coupled_field, + local_index, point_interface_data, coupled_field, self_field); specfem::point::boundary point_boundary; - specfem::assembly::load_on_device(index.self_index, boundaries, + specfem::assembly::load_on_device(index, boundaries, point_boundary); if constexpr (boundary_tag == specfem::element::boundary_tag:: acoustic_free_surface) { specfem::boundary_conditions::apply_boundary_conditions( point_boundary, self_field); } - specfem::assembly::atomic_add_on_device(index.self_index, field, + specfem::assembly::atomic_add_on_device(index, field, self_field); });