From 087ad9844b910d9b61a1a9052207c86798aa339f Mon Sep 17 00:00:00 2001 From: Lucas Sawade Date: Tue, 10 Dec 2024 16:10:38 -0500 Subject: [PATCH 1/3] Added the anisotropic mediem header --- include/medium/compute_wavefield.hpp | 1 + .../elastic_anisotropic2d.hpp | 144 ++++++++++++++++++ 2 files changed, 145 insertions(+) create mode 100644 include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp diff --git a/include/medium/compute_wavefield.hpp b/include/medium/compute_wavefield.hpp index d70616a3..9710dbf0 100644 --- a/include/medium/compute_wavefield.hpp +++ b/include/medium/compute_wavefield.hpp @@ -1,6 +1,7 @@ #pragma once #include "acoustic_isotropic2d/acoustic_isotropic2d.hpp" +#include "elastic_anisotropic2d/elastic_anisotropic2d.hpp" #include "elastic_isotropic2d/elastic_isotropic2d.hpp" #include diff --git a/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp b/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp new file mode 100644 index 00000000..3993c7ba --- /dev/null +++ b/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp @@ -0,0 +1,144 @@ +#pragma once + +#include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" +#include "enumerations/wavefield.hpp" +#include "point/field_derivatives.hpp" +#include "point/properties.hpp" +#include "point/stress.hpp" +#include + +namespace specfem { +namespace medium { + +template +KOKKOS_INLINE_FUNCTION specfem::point::stress< + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, + UseSIMD> +impl_compute_stress( + const specfem::point::properties< + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, + specfem::element::property_tag::anisotropic, UseSIMD> &properties, + const specfem::point::field_derivatives< + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, + UseSIMD> &field_derivatives) { + + using datatype = + typename specfem::datatype::simd::datatype; + const auto &du = field_derivatives.du; + + datatype sigma_xx, sigma_zz, sigma_xz; + + // P_SV case + // sigma_xx + sigma_xx = properties.c11 * du(0, 0) + properties.c13 * du(1, 1) + + properties.c15 * (du(1, 0) + du(0, 1)); + + // sigma_zz + sigma_zz = properties.c13 * du(0, 0) + properties.c33 * du(1, 1) + + properties.c35 * (du(1, 0) + du(0, 1)); + + // sigma_xz + sigma_xz = properties.c15 * du(0, 0) + properties.c35 * du(1, 1) + + properties.c55 * (du(1, 0) + du(0, 1)); + + specfem::datatype::VectorPointViewType T; + + T(0, 0) = sigma_xx; + T(0, 1) = sigma_xz; + T(1, 0) = sigma_xz; + T(1, 1) = sigma_zz; + + return { T }; +} + +template +KOKKOS_FUNCTION void impl_compute_wavefield( + const std::integral_constant, + const std::integral_constant, + const std::integral_constant, + const MemberType &team, const IteratorType &iterator, + const specfem::compute::assembly &assembly, + const QuadratureType &quadrature, const ChunkFieldType &field, + const specfem::wavefield::type wavefield_component, + WavefieldViewType wavefield) { + + using FieldDerivativesType = + specfem::point::field_derivatives; + + using PointPropertyType = specfem::point::properties< + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, + specfem::element::property_tag::isotropic, false>; + + const auto &properties = assembly.properties; + + const auto &active_field = [&]() { + if (wavefield_component == specfem::wavefield::type::displacement) { + return field.displacement; + } else if (wavefield_component == specfem::wavefield::type::velocity) { + return field.velocity; + } else if (wavefield_component == specfem::wavefield::type::acceleration) { + return field.acceleration; + } else if (wavefield_component == specfem::wavefield::type::pressure) { + return field.displacement; + } else { + Kokkos::abort("component not supported"); + } + }(); + + if (wavefield_component == specfem::wavefield::type::pressure) { + + specfem::algorithms::gradient( + team, iterator, assembly.partial_derivatives, quadrature.hprime_gll, + active_field, + [&](const typename IteratorType::index_type &iterator_index, + const FieldDerivativesType::ViewType &du) { + const auto &index = iterator_index.index; + PointPropertyType point_property; + + specfem::compute::load_on_device(index, properties, point_property); + + // P_SV case + // sigma_xx + const auto sigma_xx = point_property.c11 * du(0, 0) + + point_property.c13 * du(1, 1) + + point_property.c15 * (du(1, 0) + du(0, 1)); + + // sigma_zz + const auto sigma_zz = point_property.c13 * du(0, 0) + + point_property.c33 * du(1, 1) + + point_property.c35 * (du(1, 0) + du(0, 1)); + + // sigma_yy + const auto sigma_yy = point_property.c12 * du(0, 0) + + point_property.c23 * du(1, 1) + + point_property.c25 * (du(1, 0) + du(0, 1)); + + wavefield(index.ispec, index.iz, index.ix, 0) = + -1.0 * (sigma_xx + sigma_zz + sigma_yy) / 3.0; + }); + + return; + } + + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team, iterator.chunk_size()), [&](const int &i) { + const auto iterator_index = iterator(i); + const auto &index = iterator_index.index; + wavefield(index.ispec, index.iz, index.ix, 0) = + active_field(iterator_index.ielement, index.iz, index.ix, 0); + wavefield(index.ispec, index.iz, index.ix, 1) = + active_field(iterator_index.ielement, index.iz, index.ix, 1); + }); + + return; +} + +} // namespace medium +} // namespace specfem From adb678608ea2617f1d4fcefef629507914de3dd2 Mon Sep 17 00:00:00 2001 From: Lucas Sawade Date: Wed, 11 Dec 2024 18:16:52 -0500 Subject: [PATCH 2/3] Throwing an error for the wavefield computation. --- .../elastic_anisotropic2d.hpp | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp b/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp index 3993c7ba..94184491 100644 --- a/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp +++ b/include/medium/elastic_anisotropic2d/elastic_anisotropic2d.hpp @@ -74,7 +74,7 @@ KOKKOS_FUNCTION void impl_compute_wavefield( using PointPropertyType = specfem::point::properties< specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, - specfem::element::property_tag::isotropic, false>; + specfem::element::property_tag::anisotropic, false>; const auto &properties = assembly.properties; @@ -93,6 +93,25 @@ KOKKOS_FUNCTION void impl_compute_wavefield( }(); if (wavefield_component == specfem::wavefield::type::pressure) { + + + // over all elements and GLL points + Kokkos::parallel_for( + Kokkos::TeamThreadRange(team, iterator.chunk_size()), [&](const int &i) { + const auto iterator_index = iterator(i); + const auto &index = iterator_index.index; + + // Load the properties + PointPropertyType point_property; + specfem::compute::load_on_device(index, properties, point_property); + + // cannot compute pressure for an anisotropic material if c12 or c23 are zero + if (point_property.c12 < 1.e-7 || point_property.c23 < 1.e-7) { + Kokkos::abort("C_12 or C_23 are zero, cannot compute pressure. Check your material properties. Or, deactivate the pressure computation."); + } + + }); + specfem::algorithms::gradient( team, iterator, assembly.partial_derivatives, quadrature.hprime_gll, From 7eeb5072ce4e2643ae65002624901eb30b4693c4 Mon Sep 17 00:00:00 2001 From: Rohit Kakodkar Date: Thu, 12 Dec 2024 12:53:33 -0500 Subject: [PATCH 3/3] 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)