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 85beec248f..f43b66be4b 100644 --- a/core/specfem/assembly/assembly/dim3/assembly.cpp +++ b/core/specfem/assembly/assembly/dim3/assembly.cpp @@ -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 }; this->mpi_interfaces = { mesh.adjacency_graph, 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 f9cf518e56..e7b0a99ed3 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:: + nonconforming) { + 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/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.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/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..c25e708728 --- /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< + AccessorType::interface_tag, AccessorType::boundary_tag, + AccessorType::connection_tag, AccessorType::flux_scheme_tag>() + .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< + AccessorType::interface_tag, AccessorType::boundary_tag, + AccessorType::connection_tag, AccessorType::flux_scheme_tag>() + .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 new file mode 100644 index 0000000000..bf8fc5519e --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.hpp @@ -0,0 +1,196 @@ +#pragma once + +#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" +#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::jacobian_matrix &jacobian_matrix, + 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_nonconforming_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..a077b2ede7 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/impl/interface_container.tpp @@ -0,0 +1,163 @@ +#pragma once + +#include "interface_container.hpp" +#include "specfem/algorithms/locate_point.hpp" +#include "specfem/point/global_coordinates.hpp" +#include + +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::jacobian_matrix + &jacobian_matrix, + 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 &num_faces = self_faces.N; + const auto &weights = mesh.h_weights; + + 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); + + 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 (ipoint_i = 0; ipoint_i < nglli; ipoint_i++) { + for (ipoint_j = 0; ipoint_j < ngllj; ipoint_j++) { + specfem::point::global_coordinates global_coord( + 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( + 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 new file mode 100644 index 0000000000..d0d41372a7 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.cpp @@ -0,0 +1,23 @@ + +#include "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 &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, 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 new file mode 100644 index 0000000000..23328e2f39 --- /dev/null +++ b/core/specfem/assembly/nonconforming_interfaces/dim3/nonconforming_interfaces.hpp @@ -0,0 +1,89 @@ +#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::face, + 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 &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 = {}); + + 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 diff --git a/core/specfem/compute/impl/compute_coupling.tpp b/core/specfem/compute/impl/compute_coupling.tpp index 38d8c634db..60a78b6e3f 100644 --- a/core/specfem/compute/impl/compute_coupling.tpp +++ b/core/specfem/compute/impl/compute_coupling.tpp @@ -2,12 +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" @@ -60,9 +66,12 @@ void compute_coupling_core_weakly_conforming( const auto num_points = assembly.mesh.element_grid.ngllx; - using parallel_config = + using parallel_config = std::conditional_t< + 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; @@ -130,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); @@ -141,9 +152,12 @@ void compute_coupling_core_nonconforming( const auto num_points = assembly.mesh.element_grid.ngllx; - using parallel_config = + using parallel_config = std::conditional_t< + 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. @@ -154,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); @@ -203,40 +240,91 @@ 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); - }); + 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( + 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, 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( + local_index, point_interface_data, coupled_field, + self_field); + + specfem::point::boundary + point_boundary; + 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, 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/tag_dispatch/is_valid.hpp b/core/specfem/tag_dispatch/is_valid.hpp index 295730506d..aa8e629f49 100644 --- a/core/specfem/tag_dispatch/is_valid.hpp +++ b/core/specfem/tag_dispatch/is_valid.hpp @@ -57,7 +57,7 @@ struct TagValueTupleBase, 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)