From 7eeb5072ce4e2643ae65002624901eb30b4693c4 Mon Sep 17 00:00:00 2001 From: Rohit Kakodkar Date: Thu, 12 Dec 2024 12:53:33 -0500 Subject: [PATCH] Added template specializations for acoustic and elastic kernels class --- include/domain/impl/kernels.hpp | 176 +++++++++++++++++++++++++++- include/domain/impl/kernels.tpp | 105 ++++++++++++++--- src/domain/impl/elements/kernel.cpp | 5 + 3 files changed, 267 insertions(+), 19 deletions(-) diff --git a/include/domain/impl/kernels.hpp b/include/domain/impl/kernels.hpp index 47cd08b1..e779d1ad 100644 --- a/include/domain/impl/kernels.hpp +++ b/include/domain/impl/kernels.hpp @@ -14,8 +14,174 @@ namespace kernels { template -class kernels { +class kernels; +template +class kernels { +public: + using quadrature_point_type = qp_type; + using dimension = specfem::dimension::dimension; + kernels() = default; + + kernels(const type_real dt, const specfem::compute::assembly &assembly, + const quadrature_point_type &quadrature_points); + + /** + * @brief Compute the interaction of stiffness matrix with wavefield at a time + * step + * + * @param istep Time step + */ + inline void compute_stiffness_interaction(const int istep) const { + isotropic_elements.compute_stiffness_interaction(istep); + anisotropic_elements.compute_stiffness_interaction(istep); + isotropic_elements_dirichlet.compute_stiffness_interaction(istep); + anisotropic_elements_dirichlet.compute_stiffness_interaction(istep); + isotropic_elements_stacey.compute_stiffness_interaction(istep); + anisotropic_elements_stacey.compute_stiffness_interaction(istep); + isotropic_elements_stacey_dirichlet.compute_stiffness_interaction(istep); + anisotropic_elements_stacey_dirichlet.compute_stiffness_interaction(istep); + return; + } + + /** + * @brief Compute the mass matrix + * + * @param dt Time step + */ + inline void compute_mass_matrix(const type_real dt) const { + isotropic_elements.compute_mass_matrix(dt); + anisotropic_elements.compute_mass_matrix(dt); + isotropic_elements_dirichlet.compute_mass_matrix(dt); + anisotropic_elements_dirichlet.compute_mass_matrix(dt); + isotropic_elements_stacey.compute_mass_matrix(dt); + anisotropic_elements_stacey.compute_mass_matrix(dt); + isotropic_elements_stacey_dirichlet.compute_mass_matrix(dt); + anisotropic_elements_stacey_dirichlet.compute_mass_matrix(dt); + return; + } + + /** + * @brief Compute the interaction of source with wavefield at a time step + * + * @param timestep Time step + */ + inline void compute_source_interaction(const int timestep) const { + isotropic_sources.compute_source_interaction(timestep); + anisotropic_sources.compute_source_interaction(timestep); + return; + } + + /** + * @brief Compute seismograms at a time step + * + * @param isig_step Seismogram time step. Same if seismogram is computed at + * every time step. + */ + inline void compute_seismograms(const int &isig_step) const { + isotropic_receivers.compute_seismograms(isig_step); + anisotropic_receivers.compute_seismograms(isig_step); + return; + } + +private: + constexpr static specfem::element::boundary_tag dirichlet = + specfem::element::boundary_tag::acoustic_free_surface; + constexpr static specfem::element::boundary_tag stacey = + specfem::element::boundary_tag::stacey; + constexpr static specfem::element::boundary_tag none = + specfem::element::boundary_tag::none; + constexpr static specfem::element::boundary_tag composite_stacey_dirichlet = + specfem::element::boundary_tag::composite_stacey_dirichlet; + constexpr static specfem::element::property_tag isotropic = + specfem::element::property_tag::isotropic; + constexpr static specfem::element::property_tag anisotropic = + specfem::element::property_tag::anisotropic; + constexpr static specfem::element::medium_tag elastic = + specfem::element::medium_tag::elastic; + + constexpr static int NGLL = quadrature_point_type::NGLL; + + template + using element_kernel = specfem::domain::impl::kernels::element_kernel< + WavefieldType, DimensionType, elastic, property, boundary, + NGLL>; ///< Underlying element kernel data structure + + template + using source_kernel = specfem::domain::impl::kernels::source_kernel< + WavefieldType, DimensionType, elastic, property, + quadrature_point_type>; ///< Underlying source kernel data structure + + template + using receiver_kernel = specfem::domain::impl::kernels::receiver_kernel< + WavefieldType, DimensionType, elastic, property, + quadrature_point_type>; ///< Underlying receiver kernel data structure + + element_kernel + isotropic_elements; ///< Stiffness kernels for isotropic elements + + element_kernel + isotropic_elements_dirichlet; ///< Stiffness kernels for isotropic + ///< elements with Dirichlet boundary + ///< conditions + + element_kernel + isotropic_elements_stacey; ///< Stiffness kernels for isotropic elements + ///< with Stacey boundary conditions + + element_kernel + isotropic_elements_stacey_dirichlet; ///< Stiffness kernels for isotropic + ///< elements with Stacey and + ///< Dirichlet boundary conditions on + ///< the same element + + element_kernel + anisotropic_elements; ///< Stiffness kernels for anisotropic elements + + element_kernel + anisotropic_elements_dirichlet; ///< Stiffness kernels for anisotropic + ///< elements with Dirichlet boundary + ///< conditions + + element_kernel + anisotropic_elements_stacey; ///< Stiffness kernels for anisotropic + ///< elements with Stacey boundary conditions + + element_kernel + anisotropic_elements_stacey_dirichlet; ///< Stiffness kernels for + ///< anisotropic elements with + ///< Stacey and Dirichlet boundary + ///< conditions on the same element + + source_kernel isotropic_sources; ///< Source kernels + ///< for isotropic + ///< elements + + source_kernel + anisotropic_sources; ///< Source + ///< kernels for + ///< anisotropic + ///< elements + + receiver_kernel + isotropic_receivers; ///< Kernels for computing seismograms within + ///< isotropic elements + + receiver_kernel + anisotropic_receivers; ///< Kernels for computing seismograms within + ///< anisotropic elements +}; + +template +class kernels { public: using quadrature_point_type = qp_type; using dimension = specfem::dimension::dimension; @@ -83,6 +249,8 @@ class kernels { specfem::element::boundary_tag::composite_stacey_dirichlet; constexpr static specfem::element::property_tag isotropic = specfem::element::property_tag::isotropic; + constexpr static specfem::element::medium_tag acoustic = + specfem::element::medium_tag::acoustic; constexpr static int NGLL = quadrature_point_type::NGLL; @@ -90,19 +258,19 @@ class kernels { specfem::element::property_tag property, specfem::element::boundary_tag boundary> using element_kernel = specfem::domain::impl::kernels::element_kernel< - WavefieldType, DimensionType, medium, property, boundary, + WavefieldType, DimensionType, acoustic, property, boundary, NGLL>; ///< Underlying element kernel data structure template using source_kernel = specfem::domain::impl::kernels::source_kernel< - WavefieldType, DimensionType, medium, property, + WavefieldType, DimensionType, acoustic, property, quadrature_point_type>; ///< Underlying source kernel data structure template using receiver_kernel = specfem::domain::impl::kernels::receiver_kernel< - WavefieldType, DimensionType, medium, property, + WavefieldType, DimensionType, acoustic, property, quadrature_point_type>; ///< Underlying receiver kernel data structure element_kernel diff --git a/include/domain/impl/kernels.tpp b/include/domain/impl/kernels.tpp index e07257e6..b8ee0c4f 100644 --- a/include/domain/impl/kernels.tpp +++ b/include/domain/impl/kernels.tpp @@ -73,15 +73,16 @@ void allocate_elements( } } - if constexpr (wavefield_type == specfem::wavefield::simulation_field::forward || - wavefield_type == specfem::wavefield::simulation_field::adjoint) { + if constexpr (wavefield_type == + specfem::wavefield::simulation_field::forward || + wavefield_type == + specfem::wavefield::simulation_field::adjoint) { std::cout << " - Element type: \n" << " - dimension : " << dimension::to_string() << "\n" - << " - Element type : " << specfem::element::to_string( - medium_tag, property_tag) - << "\n" + << " - Element type : " + << specfem::element::to_string(medium_tag, property_tag) << "\n" << " - Boundary Conditions : " << specfem::domain::impl::boundary_conditions::print_boundary_tag< boundary_tag>() @@ -190,12 +191,85 @@ void allocate_isotropic_receivers( } // namespace template -specfem::domain::impl::kernels:: - kernels::kernels( - const type_real dt, const specfem::compute::assembly &assembly, - const qp_type &quadrature_points) { + specfem::dimension::type DimensionType, typename qp_type> +specfem::domain::impl::kernels::kernels< + WavefieldType, DimensionType, specfem::element::medium_tag::elastic, + qp_type>::kernels(const type_real dt, + const specfem::compute::assembly &assembly, + const qp_type &quadrature_points) { + + const int nspec = assembly.mesh.nspec; + specfem::kokkos::HostView1d element_tags( + "specfem::domain::domain::element_tag", nspec); + + // ----------------------------------------------------------- + for (int ispec = 0; ispec < nspec; ispec++) { + element_tags(ispec) = + element_tag(assembly.properties.h_element_types(ispec), + assembly.properties.h_element_property(ispec), + assembly.boundaries.boundary_tags(ispec)); + } + + if constexpr (WavefieldType == + specfem::wavefield::simulation_field::forward || + WavefieldType == + specfem::wavefield::simulation_field::adjoint) { + std::cout << " Element Statistics \n" + << "------------------------------\n" + << "- Types of elements in elastic medium :\n\n"; + } + + // ----------------------------------------------------------- + + // Allocate isotropic elements with dirichlet boundary conditions + allocate_elements(assembly, element_tags, isotropic_elements_dirichlet); + + // Allocate aniostropic elements with dirichlet boundary conditions + allocate_elements(assembly, element_tags, anisotropic_elements_dirichlet); + + // Allocate isotropic elements with stacey boundary conditions + allocate_elements(assembly, element_tags, isotropic_elements_stacey); + + // Allocate anisotropic elements with stacey boundary conditions + allocate_elements(assembly, element_tags, anisotropic_elements_stacey); + + // Allocate isotropic elements with stacey dirichlet boundary conditions + allocate_elements(assembly, element_tags, + isotropic_elements_stacey_dirichlet); + + // Allocate anisotropic elements with stacey dirichlet boundary conditions + allocate_elements(assembly, element_tags, + anisotropic_elements_stacey_dirichlet); + + // Allocate isotropic elements + allocate_elements(assembly, element_tags, isotropic_elements); + + // Allocate anisotropic elements + allocate_elements(assembly, element_tags, anisotropic_elements); + + // Allocate isotropic sources + + allocate_isotropic_sources(assembly, quadrature_points, isotropic_sources); + + // Allocate isotropic receivers + + allocate_isotropic_receivers(assembly, quadrature_points, + isotropic_receivers); + + // Compute mass matrices + + this->compute_mass_matrix(dt); + + return; +} + +template +specfem::domain::impl::kernels::kernels< + WavefieldType, DimensionType, specfem::element::medium_tag::acoustic, + qp_type>::kernels(const type_real dt, + const specfem::compute::assembly &assembly, + const qp_type &quadrature_points) { const int nspec = assembly.mesh.nspec; specfem::kokkos::HostView1d element_tags( @@ -209,12 +283,13 @@ specfem::domain::impl::kernels:: assembly.boundaries.boundary_tags(ispec)); } - if constexpr (WavefieldType == specfem::wavefield::simulation_field::forward || - WavefieldType == specfem::wavefield::simulation_field::adjoint) { + if constexpr (WavefieldType == + specfem::wavefield::simulation_field::forward || + WavefieldType == + specfem::wavefield::simulation_field::adjoint) { std::cout << " Element Statistics \n" << "------------------------------\n" - << "- Types of elements in " << specfem::element::to_string(medium) - << " medium :\n\n"; + << "- Types of elements in acoustic medium :\n\n"; } // ----------------------------------------------------------- diff --git a/src/domain/impl/elements/kernel.cpp b/src/domain/impl/elements/kernel.cpp index bfe2d824..9b6435b9 100644 --- a/src/domain/impl/elements/kernel.cpp +++ b/src/domain/impl/elements/kernel.cpp @@ -11,6 +11,7 @@ constexpr static auto elastic = specfem::element::medium_tag::elastic; constexpr static auto acoustic = specfem::element::medium_tag::acoustic; constexpr static auto isotropic = specfem::element::property_tag::isotropic; +constexpr static auto anisotropic = specfem::element::property_tag::anisotropic; constexpr static auto dirichlet = specfem::element::boundary_tag::acoustic_free_surface; @@ -79,8 +80,12 @@ constexpr static auto composite_stacey_dirichlet = GENERATE_KERNELS(elastic, isotropic, 5) +GENERATE_KERNELS(elastic, anisotropic, 5) + GENERATE_KERNELS(acoustic, isotropic, 5) GENERATE_KERNELS(elastic, isotropic, 8) +GENERATE_KERNELS(elastic, anisotropic, 8) + GENERATE_KERNELS(acoustic, isotropic, 8)