diff --git a/core/specfem/assembly/assembly.hpp b/core/specfem/assembly/assembly.hpp index 0b7325938..ad8a48d1a 100644 --- a/core/specfem/assembly/assembly.hpp +++ b/core/specfem/assembly/assembly.hpp @@ -16,7 +16,7 @@ * specfem::assembly::mesh, @ref specfem::assembly::jacobian_matrix, etc) that * store data computed at all GLL points * - Data access functions : Each container provides functions to load/store - * data on device/host (e.g., @c load_on_device , @c store_on_device ) + * data on device/host (e.g., @c lod_on_device , @c store_on_device ) * */ namespace specfem::assembly { diff --git a/core/specfem/assembly/assembly/dim2/assembly.cpp b/core/specfem/assembly/assembly/dim2/assembly.cpp index 12c741221..831042162 100644 --- a/core/specfem/assembly/assembly/dim2/assembly.cpp +++ b/core/specfem/assembly/assembly/dim2/assembly.cpp @@ -29,6 +29,13 @@ specfem::assembly::assembly::assembly( this->mesh.element_grid.ngllz, this->mesh, this->element_types, flux_scheme_config }; this->jacobian_matrix = { this->mesh }; + + this->field_derivative_storage = { this->element_types, this->mesh.nspec, + this->mesh.element_grid.ngllz, + this->mesh.element_grid.ngllx }; + + // this->attenuation = { reference}; + this->properties = { this->mesh.nspec, this->mesh.element_grid.ngllz, this->mesh.element_grid.ngllx, diff --git a/core/specfem/assembly/assembly/dim2/assembly.hpp b/core/specfem/assembly/assembly/dim2/assembly.hpp index 8a280c303..9e67d3e72 100644 --- a/core/specfem/assembly/assembly/dim2/assembly.hpp +++ b/core/specfem/assembly/assembly/dim2/assembly.hpp @@ -1,10 +1,12 @@ #pragma once +#include "specfem/assembly/attenuation.hpp" #include "specfem/assembly/boundaries.hpp" #include "specfem/assembly/boundary_values.hpp" #include "specfem/assembly/conforming_interfaces.hpp" #include "specfem/assembly/element_intersections.hpp" #include "specfem/assembly/element_types.hpp" +#include "specfem/assembly/field_derivative_storage.hpp" #include "specfem/assembly/fields.hpp" #include "specfem/assembly/info.hpp" #include "specfem/assembly/jacobian_matrix.hpp" @@ -78,6 +80,17 @@ template <> struct assembly { */ specfem::assembly::jacobian_matrix jacobian_matrix; + /** + * @brief Storage for field derivatives at every quadrature point for elements + * that require it (e.g., attenuating elements). + */ + specfem::assembly::FieldDerivativeStorage + field_derivative_storage; + /** + * @brief Attenuation properties for the mesh at every quadrature point + */ + // specfem::assembly::attenuation attenuation; + /** * @brief Material properties for the mesh at every quadrature point * diff --git a/core/specfem/assembly/assembly/dim3/assembly.cpp b/core/specfem/assembly/assembly/dim3/assembly.cpp index 9da587bdd..ea6abceb5 100644 --- a/core/specfem/assembly/assembly/dim3/assembly.cpp +++ b/core/specfem/assembly/assembly/dim3/assembly.cpp @@ -41,6 +41,11 @@ specfem::assembly::assembly::assembly( this->jacobian_matrix = { this->mesh }; + this->field_derivative_storage = { this->element_types, this->mesh.nspec, + this->mesh.element_grid.ngllz, + this->mesh.element_grid.nglly, + this->mesh.element_grid.ngllx }; + this->properties = { nspec, ngllz, nglly, ngllx, mesh.materials, this->element_types }; diff --git a/core/specfem/assembly/assembly/dim3/assembly.hpp b/core/specfem/assembly/assembly/dim3/assembly.hpp index 32c475de9..2702bbd2d 100644 --- a/core/specfem/assembly/assembly/dim3/assembly.hpp +++ b/core/specfem/assembly/assembly/dim3/assembly.hpp @@ -1,9 +1,11 @@ #pragma once +#include "specfem/assembly/attenuation.hpp" #include "specfem/assembly/boundaries.hpp" #include "specfem/assembly/boundary_values.hpp" #include "specfem/assembly/conforming_interfaces.hpp" #include "specfem/assembly/element_intersections.hpp" +#include "specfem/assembly/field_derivative_storage.hpp" #include "specfem/assembly/fields.hpp" #include "specfem/assembly/info.hpp" #include "specfem/assembly/jacobian_matrix.hpp" @@ -71,6 +73,18 @@ template <> struct assembly { */ specfem::assembly::jacobian_matrix jacobian_matrix; + /** + * @brief Storage for field derivatives at every quadrature point for elements + * that require it (e.g., attenuating elements). + */ + specfem::assembly::FieldDerivativeStorage + field_derivative_storage; + + /** + * @brief Attenuation properties for the mesh at every quadrature point + */ + // specfem::assembly::attenuation attenuation; + /** * @brief Material properties for the mesh at every quadrature point * diff --git a/core/specfem/assembly/attenuation.hpp b/core/specfem/assembly/attenuation.hpp new file mode 100644 index 000000000..aaaa2ebe0 --- /dev/null +++ b/core/specfem/assembly/attenuation.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include "specfem/enums.hpp" + +namespace specfem::assembly { + +template struct Attenuation; + +} // namespace specfem::assembly + +#include "specfem/assembly/attenuation/dim2/attenuation.hpp" +#include "specfem/assembly/attenuation/dim3/attenuation.hpp" +#include "specfem/assembly/attenuation/load_on_device.hpp" +#include "specfem/assembly/attenuation/store_on_device.hpp" diff --git a/core/specfem/assembly/attenuation/dim2/attenuation.cpp b/core/specfem/assembly/attenuation/dim2/attenuation.cpp new file mode 100644 index 000000000..4b041b5f7 --- /dev/null +++ b/core/specfem/assembly/attenuation/dim2/attenuation.cpp @@ -0,0 +1,87 @@ +#include "specfem/assembly/attenuation/dim2/attenuation.hpp" +#include "specfem/attenuation.hpp" +#include "specfem/setup.hpp" +#include "specfem/utilities/band.hpp" + +#include + +void specfem::assembly::Attenuation:: + init_memory_variables( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim2> &element_types, + const specfem::assembly::mesh + &mesh, + const specfem::mesh::materials + &materials, + const specfem::units::Hertz fc, const specfem::units::Hertz f0, + const specfem::utilities::Band &band, + const Kokkos::View + &tau_sigma){ + + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV), + PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { + auto elements = element_types.get_elements_on_host( + _medium_tag_, _property_tag_, _attenuation_tag_); + _attn_medium_ = { elements, mesh, materials, ngllz, ngllx, + fc, f0, band, tau_sigma }; + }) + } + +specfem::assembly::Attenuation:: + Attenuation( + const specfem::units::Hertz reference_frequency, + const specfem::utilities::Band band_in, + const bool auto_compute_attenuation_band, const type_real deltat, + const specfem::assembly::mesh + &mesh, + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim2> &element_types, + const specfem::assembly::Info + &info, + const specfem::mesh::materials + &materials) + : ngllz(mesh.element_grid.ngllz), ngllx(mesh.element_grid.ngllx), + nspec(mesh.nspec), f0(reference_frequency), + auto_compute_attenuation_band(auto_compute_attenuation_band) { + + specfem::utilities::Band band = + auto_compute_attenuation_band + ? attenuation::compute_band( + specfem::units::Seconds(info.largest_minimum_period)) + : band_in; + + using specfem::units::unit_symbols::Hz; + + // Compute the band center frequency for attenuation scaling (logarithmic + // center) + auto fc = + specfem::utilities::logarithmic_center(band.min.raw(), band.max.raw()) * + Hz; + + // Compute tau_sigma once (shared across all elements) + auto tau_sigma = specfem::attenuation::compute_tau_sigma(band); + + // Compute Runge-Kutta memory-variable update coefficients + auto rk = specfem::attenuation::compute_integration_factors(tau_sigma, + deltat); + + this->alpha_rk = rk.alpha; + this->beta_rk = rk.beta; + this->gamma_rk = rk.gamma; + + // Copy RK coefficients to device for use in device kernels + this->d_alpha_rk = Kokkos::View("d_alpha_rk"); + Kokkos::deep_copy(this->d_alpha_rk, this->alpha_rk); + this->d_beta_rk = Kokkos::View("d_beta_rk"); + Kokkos::deep_copy(this->d_beta_rk, this->beta_rk); + this->d_gamma_rk = Kokkos::View("d_gamma_rk"); + Kokkos::deep_copy(this->d_gamma_rk, this->gamma_rk); + + init_memory_variables(element_types, mesh, materials, fc, f0, band, + tau_sigma); +} diff --git a/core/specfem/assembly/attenuation/dim2/attenuation.hpp b/core/specfem/assembly/attenuation/dim2/attenuation.hpp new file mode 100644 index 000000000..8059212c6 --- /dev/null +++ b/core/specfem/assembly/attenuation/dim2/attenuation.hpp @@ -0,0 +1,171 @@ +#pragma once + +#include "specfem/assembly/attenuation.hpp" +#include "specfem/assembly/attenuation/impl/attenuation_medium.hpp" +#include "specfem/assembly/element_types.hpp" +#include "specfem/assembly/info.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/constants.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros.hpp" +#include "specfem/mesh/dim2/materials/materials.hpp" +#include "specfem/setup.hpp" +#include + +namespace specfem::assembly { + +template <> +struct Attenuation + : public specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim2> { + + /** + * @name Type Definitions + * + */ + ///@{ + + /** + * @brief Base container type providing data access infrastructure + * + * @see specfem::data_access::Container + */ + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim2>; + + /** + * @brief Kokkos view type for per-element index mapping arrays + */ + using IndexViewType = Kokkos::View; + ///@} + + /** + * @name Compile Time Definitions + * + */ + ///@{ + + constexpr static auto dimension_tag = + specfem::element::dimension_tag::dim2; ///< Dimension tag for 2D + constexpr static int N_SLS = + specfem::constants::N_SLS; ///< Number of standard linear solids + + ///@} + + // Number of GLL points and elements + int ngllz; ///< Number of GLL points in the z-direction + int ngllx; ///< Number of GLL points in the x-direction + int nspec; ///< Total number of spectral elements + + // Attenuation parameters + specfem::units::Hertz f0; ///< Reference frequency + bool auto_compute_attenuation_band; ///< Whether to auto-compute the + ///< attenuation band + + // Runge-Kutta attenuation factors (one coefficient per SLS mechanism) + Kokkos::View + alpha_rk; + Kokkos::View + beta_rk; + Kokkos::View + gamma_rk; + + // Device-accessible Runge-Kutta factors for use in device kernels + Kokkos::View + d_alpha_rk; + Kokkos::View + d_beta_rk; + Kokkos::View + d_gamma_rk; + + // One attenuation_medium member per (medium, property) combination + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV), + PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + DECLARE(((specfem::assembly::impl::attenuation_medium, + (_DIMENSION_TAG_, _MEDIUM_TAG_, _PROPERTY_TAG_, + _ATTENUATION_TAG_)), + attn_medium))) + + Attenuation() = default; + + Attenuation( + const specfem::units::Hertz reference_frequency, + const specfem::utilities::Band band_in, + const bool auto_compute_attenuation_band, const type_real deltat, + const specfem::assembly::mesh + &mesh, + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim2> &element_types, + const specfem::assembly::Info + &info, + const specfem::mesh::materials + &materials); + + void init_memory_variables( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim2> &element_types, + const specfem::assembly::mesh + &mesh, + const specfem::mesh::materials + &materials, + const specfem::units::Hertz fc, const specfem::units::Hertz f0, + const specfem::utilities::Band &band, + const Kokkos::View + &tau_sigma); + + /** + * @brief Access the attenuation_medium for a given medium and property. + */ + template + KOKKOS_INLINE_FUNCTION constexpr specfem::assembly::impl::attenuation_medium< + specfem::element::dimension_tag::dim2, MediumTag, PropertyTag, + specfem::element::attenuation_tag::constant_isotropic> const & + get_container() const { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV), + PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { + if constexpr (_medium_tag_ == MediumTag && + _property_tag_ == PropertyTag) { + return _attn_medium_; + } + }) + Kokkos::abort("Invalid medium type detected in attenuation"); + SUPPRESS_UNREACHABLE(return {};) + } + + /** + * @brief Copy attenuation data to host mirrors. + */ + void copy_to_host() { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV), + PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { _attn_medium_.copy_to_host(); }) + } + + /** + * @brief Copy attenuation data to device views. + */ + void copy_to_device() { + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV), PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { _attn_medium_.copy_to_device(); }) + } +}; + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/attenuation/dim3/attenuation.cpp b/core/specfem/assembly/attenuation/dim3/attenuation.cpp new file mode 100644 index 000000000..bd13be4fd --- /dev/null +++ b/core/specfem/assembly/attenuation/dim3/attenuation.cpp @@ -0,0 +1,85 @@ +#include "specfem/assembly/attenuation/dim3/attenuation.hpp" +#include "specfem/attenuation.hpp" +#include "specfem/utilities.hpp" +#include + +void specfem::assembly::Attenuation:: + init_memory_variables( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim3> &element_types, + const specfem::assembly::mesh + &mesh, + const specfem::mesh::materials + &materials, + const specfem::units::Hertz fc, const specfem::units::Hertz f0, + const specfem::utilities::Band &band, + const Kokkos::View + &tau_sigma){ + + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { + auto elements = element_types.get_elements_on_host( + _medium_tag_, _property_tag_, _attenuation_tag_); + _attn_medium_ = { elements, mesh, materials, ngllz, nglly, + ngllx, fc, f0, band, tau_sigma }; + }) + } + +specfem::assembly::Attenuation:: + Attenuation( + const specfem::units::Hertz reference_frequency, + const specfem::utilities::Band band_in, + const bool auto_compute_attenuation_band, const type_real deltat, + const specfem::assembly::mesh + &mesh, + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim3> &element_types, + const specfem::assembly::Info + &info, + const specfem::mesh::materials + &materials) + : ngllz(mesh.element_grid.ngllz), nglly(mesh.element_grid.nglly), + ngllx(mesh.element_grid.ngllx), nspec(mesh.nspec), + f0(reference_frequency), + auto_compute_attenuation_band(auto_compute_attenuation_band) { + + specfem::utilities::Band band = + auto_compute_attenuation_band + ? attenuation::compute_band( + specfem::units::Seconds(info.largest_minimum_period)) + : band_in; + + using specfem::units::unit_symbols::Hz; + + // Compute the band center frequency for attenuation scaling (logarithmic + // center) + auto fc = + specfem::utilities::logarithmic_center(band.min.raw(), band.max.raw()) * + Hz; + + // Compute tau_sigma once (shared across all elements) + auto tau_sigma = specfem::attenuation::compute_tau_sigma(band); + + // Compute Runge-Kutta memory-variable update coefficients + auto rk = specfem::attenuation::compute_integration_factors(tau_sigma, + deltat); + this->alpha_rk = rk.alpha; + this->beta_rk = rk.beta; + this->gamma_rk = rk.gamma; + + // Copy RK coefficients to device for use in device kernels + this->d_alpha_rk = Kokkos::View("d_alpha_rk"); + Kokkos::deep_copy(this->d_alpha_rk, this->alpha_rk); + this->d_beta_rk = Kokkos::View("d_beta_rk"); + Kokkos::deep_copy(this->d_beta_rk, this->beta_rk); + this->d_gamma_rk = Kokkos::View("d_gamma_rk"); + Kokkos::deep_copy(this->d_gamma_rk, this->gamma_rk); + + init_memory_variables(element_types, mesh, materials, fc, f0, band, + tau_sigma); +} diff --git a/core/specfem/assembly/attenuation/dim3/attenuation.hpp b/core/specfem/assembly/attenuation/dim3/attenuation.hpp new file mode 100644 index 000000000..ef543c111 --- /dev/null +++ b/core/specfem/assembly/attenuation/dim3/attenuation.hpp @@ -0,0 +1,173 @@ +#pragma once + +#include "specfem/assembly/attenuation.hpp" +#include "specfem/assembly/attenuation/impl/attenuation_medium.hpp" +#include "specfem/assembly/element_types.hpp" +#include "specfem/assembly/info.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/constants.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros.hpp" +#include "specfem/mesh/dim3/materials/materials.hpp" +#include "specfem/setup.hpp" +#include "specfem/units.hpp" +#include + +namespace specfem::assembly { + +template <> +struct Attenuation + : public specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim3> { + + /** + * @name Type Definitions + * + */ + ///@{ + + /** + * @brief Base container type providing data access infrastructure + * + * @see specfem::data_access::Container + */ + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim3>; + + /** + * @brief Kokkos view type for per-element index mapping arrays + */ + using IndexViewType = Kokkos::View; + ///@} + + /** + * @name Compile Time Definitions + * + */ + ///@{ + + constexpr static auto dimension_tag = + specfem::element::dimension_tag::dim3; ///< Dimension tag for 3D + constexpr static int N_SLS = + specfem::constants::N_SLS; ///< Number of standard linear solids + + ///@} + + // Number of GLL points and elements + int ngllz; ///< Number of GLL points in the z-direction + int nglly; ///< Number of GLL points in the y-direction + int ngllx; ///< Number of GLL points in the x-direction + int nspec; ///< Total number of spectral elements + + // Attenuation parameters + specfem::units::Hertz f0; ///< Reference frequency + bool auto_compute_attenuation_band; ///< Whether to auto-compute the + ///< attenuation band + + // Runge-Kutta attenuation factors (one coefficient per SLS mechanism) + Kokkos::View + alpha_rk; + Kokkos::View + beta_rk; + Kokkos::View + gamma_rk; + + // Device-accessible Runge-Kutta factors for use in device kernels + Kokkos::View + d_alpha_rk; + Kokkos::View + d_beta_rk; + Kokkos::View + d_gamma_rk; + + // One attenuation_medium member per (medium, property) combination + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), + PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + DECLARE(((specfem::assembly::impl::attenuation_medium, + (_DIMENSION_TAG_, _MEDIUM_TAG_, _PROPERTY_TAG_, + _ATTENUATION_TAG_)), + attn_medium))) + + Attenuation() = default; + + Attenuation( + const specfem::units::Hertz reference_frequency, + const specfem::utilities::Band band_in, + const bool auto_compute_attenuation_band, const type_real deltat, + const specfem::assembly::mesh + &mesh, + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim3> &element_types, + const specfem::assembly::Info + &info, + const specfem::mesh::materials + &materials); + + void init_memory_variables( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim3> &element_types, + const specfem::assembly::mesh + &mesh, + const specfem::mesh::materials + &materials, + const specfem::units::Hertz fc, const specfem::units::Hertz f0, + const specfem::utilities::Band &band, + const Kokkos::View + &tau_sigma); + + /** + * @brief Access the attenuation_medium for a given medium and property. + */ + template + KOKKOS_INLINE_FUNCTION constexpr specfem::assembly::impl::attenuation_medium< + specfem::element::dimension_tag::dim3, MediumTag, PropertyTag, + specfem::element::attenuation_tag::constant_isotropic> const & + get_container() const { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), + PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { + if constexpr (_medium_tag_ == MediumTag && + _property_tag_ == PropertyTag) { + return _attn_medium_; + } + }) + Kokkos::abort("Invalid medium type detected in attenuation"); + SUPPRESS_UNREACHABLE(return {};) + } + + /** + * @brief Copy attenuation data to host mirrors. + */ + void copy_to_host() { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), + PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { _attn_medium_.copy_to_host(); }) + } + + /** + * @brief Copy attenuation data to device views. + */ + void copy_to_device() { + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC), PROPERTY_TAG(ISOTROPIC), + ATTENUATION_TAG(CONSTANT_ISOTROPIC)), + CAPTURE(attn_medium) { _attn_medium_.copy_to_device(); }) + } +}; + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/attenuation/impl/attenuation_medium.hpp b/core/specfem/assembly/attenuation/impl/attenuation_medium.hpp new file mode 100644 index 000000000..34cc81339 --- /dev/null +++ b/core/specfem/assembly/attenuation/impl/attenuation_medium.hpp @@ -0,0 +1,16 @@ +#pragma once + +#include "specfem/enums.hpp" + +namespace specfem::assembly::impl { + +template +struct attenuation_medium; + +} // namespace specfem::assembly::impl + +#include "specfem/assembly/attenuation/impl/attenuation_medium/dim2/constant_isotropic.tpp" +#include "specfem/assembly/attenuation/impl/attenuation_medium/dim3/constant_isotropic.tpp" diff --git a/core/specfem/assembly/attenuation/impl/attenuation_medium/dim2/constant_isotropic.tpp b/core/specfem/assembly/attenuation/impl/attenuation_medium/dim2/constant_isotropic.tpp new file mode 100644 index 000000000..bf520539b --- /dev/null +++ b/core/specfem/assembly/attenuation/impl/attenuation_medium/dim2/constant_isotropic.tpp @@ -0,0 +1,267 @@ +#pragma once + +#include "specfem/assembly/attenuation/impl/attenuation_medium.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/attenuation/compute_factors.hpp" +#include "specfem/attenuation/compute_tau_eps.hpp" +#include "specfem/attenuation/compute_tau_sigma.hpp" +#include "specfem/constants.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/enums.hpp" +#include "specfem/mesh/dim2/materials/materials.hpp" +#include "specfem/setup.hpp" +#include "specfem/utilities/logarithmic_center.hpp" +#include "specfem/units.hpp" +#include + +namespace specfem::assembly::impl { + +template +struct attenuation_medium : + specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim2> { + + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim2>; + + using view_type = typename base_type::vector_type< + type_real, Kokkos::DefaultExecutionSpace::memory_space>; + + constexpr static int N_SLS = specfem::constants::N_SLS; + + + // Host-only per-element scale factors + Kokkos::View h_kappa_scale; + Kokkos::View h_mu_scale; + + // Views: shape [nspec_attn][ngllz][ngllx][N_SLS] + view_type kappa_relaxation_rate; + view_type::HostMirror h_kappa_relaxation_rate; + view_type mu_relaxation_rate; + view_type::HostMirror h_mu_relaxation_rate; + view_type memory_variable_kappa; + view_type::HostMirror h_memory_variable_kappa; + view_type memory_variable_Rxx; + view_type::HostMirror h_memory_variable_Rxx; + view_type memory_variable_Rxz; + view_type::HostMirror h_memory_variable_Rxz; + + // Index mapping: global ispec -> compact attenuation index (-1 if not attenuating) + Kokkos::View + h_attenuation_index_mapping; + Kokkos::View + attenuation_index_mapping; + + attenuation_medium() = default; + + attenuation_medium( + const Kokkos::View elements, + const specfem::assembly::mesh + &mesh, + const specfem::mesh::materials + &materials, + const int ngllz, const int ngllx, const specfem::units::Hertz fc, const specfem::units::Hertz f0, + const specfem::utilities::Band &band, + const Kokkos::View &tau_sigma) { + + const int nspec_attn = elements.extent(0); + + // 1. Allocate all views + h_kappa_scale = + Kokkos::View( + "kappa_scale", nspec_attn); + h_mu_scale = + Kokkos::View( + "mu_scale", nspec_attn); + kappa_relaxation_rate = + view_type("kappa_relaxation_rate", nspec_attn, ngllz, ngllx, N_SLS); + h_kappa_relaxation_rate = + specfem::datatype::create_mirror_view(kappa_relaxation_rate); + mu_relaxation_rate = + view_type("mu_relaxation_rate", nspec_attn, ngllz, ngllx, N_SLS); + h_mu_relaxation_rate = + specfem::datatype::create_mirror_view(mu_relaxation_rate); + memory_variable_kappa = + view_type("mem_kappa", nspec_attn, ngllz, ngllx, N_SLS); + h_memory_variable_kappa = + specfem::datatype::create_mirror_view(memory_variable_kappa); + memory_variable_Rxx = + view_type("mem_Rxx", nspec_attn, ngllz, ngllx, N_SLS); + h_memory_variable_Rxx = + specfem::datatype::create_mirror_view(memory_variable_Rxx); + memory_variable_Rxz = + view_type("mem_Rxz", nspec_attn, ngllz, ngllx, N_SLS); + h_memory_variable_Rxz = + specfem::datatype::create_mirror_view(memory_variable_Rxz); + + // Allocate and populate the inverse index mapping (global ispec -> compact index) + h_attenuation_index_mapping = + Kokkos::View( + "h_attenuation_index_mapping", mesh.nspec); + Kokkos::deep_copy(h_attenuation_index_mapping, -1); + for (int i = 0; i < nspec_attn; ++i) { + h_attenuation_index_mapping(elements(i)) = i; + } + attenuation_index_mapping = + Kokkos::View( + "attenuation_index_mapping", mesh.nspec); + Kokkos::deep_copy(attenuation_index_mapping, h_attenuation_index_mapping); + + // Sync zero-initialized device views to host mirrors + copy_to_host(); + + if (nspec_attn == 0) { + return; + } + + // Sanity check input frequencies before looping over elements + if (fc.raw() <= 0 || f0.raw() <= 0) { + throw std::runtime_error( + "Center frequency fc and reference frequency f0 must be positive."); + } + + // 3. Loop over elements + for (int i = 0; i < nspec_attn; ++i) { + const int ispec = elements(i); + const int mesh_ispec = mesh.compute_to_mesh(ispec); + + auto material = materials.template get_material< + specfem::element::medium_tag::elastic_psv, PropertyTag, + specfem::element::attenuation_tag::constant_isotropic>(mesh_ispec); + + auto computed_values = material.compute_attenuation_properties(fc.raw(), f0.raw(), band, tau_sigma); + + auto kappa_props = computed_values.kappa_attenuation_properties; + auto mu_props = computed_values.mu_attenuation_properties; + + // Per-GLL fill + for (int j = 0; j < N_SLS; ++j) { + + // Compute per element relaxation rates + const type_real tauinv_j = 1.0 / tau_sigma(j); + auto kappa_rr_j = kappa_props.beta(j) * tauinv_j / kappa_props.one_minus_sum_beta; + auto mu_rr_j = 2.0 * mu_props.beta(j) * tauinv_j / mu_props.one_minus_sum_beta; + + // Assigning relaxation rates to all GLL points + for (int iz = 0; iz < ngllz; ++iz) { + for (int ix = 0; ix < ngllx; ++ix) { + h_kappa_relaxation_rate(i, iz, ix, j) = kappa_rr_j; + h_mu_relaxation_rate(i, iz, ix, j) = mu_rr_j; + } + } + } + } + + + // 4. Push all host data (kappa/mu_cf filled; memory variables zero) to device + copy_to_device(); + } + + void copy_to_host() { + Kokkos::deep_copy(h_kappa_relaxation_rate, kappa_relaxation_rate); + Kokkos::deep_copy(h_mu_relaxation_rate, mu_relaxation_rate); + Kokkos::deep_copy(h_memory_variable_kappa, memory_variable_kappa); + Kokkos::deep_copy(h_memory_variable_Rxx, memory_variable_Rxx); + Kokkos::deep_copy(h_memory_variable_Rxz, memory_variable_Rxz); + } + + void copy_to_device() { + Kokkos::deep_copy(kappa_relaxation_rate, h_kappa_relaxation_rate); + Kokkos::deep_copy(mu_relaxation_rate, h_mu_relaxation_rate); + Kokkos::deep_copy(memory_variable_kappa, h_memory_variable_kappa); + Kokkos::deep_copy(memory_variable_Rxx, h_memory_variable_Rxx); + Kokkos::deep_copy(memory_variable_Rxz, h_memory_variable_Rxz); + } + + /** + * @brief Load attenuation data for a single GLL point from device views + * into a point-local attenuation struct. + * + * Populates relaxation rates and memory variables. The global RK + * coefficients are NOT populated here; they are added by the outer + * load_on_device free function. + */ + template + KOKKOS_INLINE_FUNCTION void load_device_values(const IndexType &index, + PointType &point) const { + const int i = attenuation_index_mapping(index.ispec); + if constexpr (!IndexType::using_simd) { + for (int j = 0; j < N_SLS; ++j) { + point.kappa_relaxation_rate(j) = + kappa_relaxation_rate(i, index.iz, index.ix, j); + point.mu_relaxation_rate(j) = + mu_relaxation_rate(i, index.iz, index.ix, j); + point.Rxx(j) = memory_variable_Rxx(i, index.iz, index.ix, j); + point.Rxz(j) = memory_variable_Rxz(i, index.iz, index.ix, j); + point.Rkappa(j) = memory_variable_kappa(i, index.iz, index.ix, j); + } + } else { + using simd = typename PointType::simd; + using mask_type = typename simd::mask_type; + using tag_type = typename simd::tag_type; + const auto mask = index.template get_mask(); + for (int j = 0; j < N_SLS; ++j) { + Kokkos::Experimental::where(mask, point.kappa_relaxation_rate(j)) + .copy_from(&kappa_relaxation_rate(i, index.iz, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.mu_relaxation_rate(j)) + .copy_from(&mu_relaxation_rate(i, index.iz, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxx(j)) + .copy_from(&memory_variable_Rxx(i, index.iz, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxz(j)) + .copy_from(&memory_variable_Rxz(i, index.iz, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rkappa(j)) + .copy_from(&memory_variable_kappa(i, index.iz, index.ix, j), + tag_type()); + } + } + } + + /** + * @brief Store evolved SLS memory variables from a point-local struct back + * to the device views. + * + * Only the memory variables (Rxx, Rxz, Rkappa) are written; relaxation + * rates are simulation-lifetime constants and are not written back. + */ + template + KOKKOS_INLINE_FUNCTION void + store_device_values(const IndexType &index, const PointType &point) const { + const int i = attenuation_index_mapping(index.ispec); + if constexpr (!IndexType::using_simd) { + for (int j = 0; j < N_SLS; ++j) { + memory_variable_Rxx(i, index.iz, index.ix, j) = point.Rxx(j); + memory_variable_Rxz(i, index.iz, index.ix, j) = point.Rxz(j); + memory_variable_kappa(i, index.iz, index.ix, j) = point.Rkappa(j); + } + } else { + using simd = typename PointType::simd; + using mask_type = typename simd::mask_type; + using tag_type = typename simd::tag_type; + const auto mask = index.template get_mask(); + for (int j = 0; j < N_SLS; ++j) { + Kokkos::Experimental::where(mask, point.Rxx(j)) + .copy_to(&memory_variable_Rxx(i, index.iz, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxz(j)) + .copy_to(&memory_variable_Rxz(i, index.iz, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rkappa(j)) + .copy_to(&memory_variable_kappa(i, index.iz, index.ix, j), + tag_type()); + } + } + } +}; + +} // namespace specfem::assembly::impl diff --git a/core/specfem/assembly/attenuation/impl/attenuation_medium/dim3/constant_isotropic.tpp b/core/specfem/assembly/attenuation/impl/attenuation_medium/dim3/constant_isotropic.tpp new file mode 100644 index 000000000..d700ece74 --- /dev/null +++ b/core/specfem/assembly/attenuation/impl/attenuation_medium/dim3/constant_isotropic.tpp @@ -0,0 +1,333 @@ +#pragma once + +#include "specfem/assembly/attenuation/impl/attenuation_medium.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/attenuation/compute_factors.hpp" +#include "specfem/attenuation/compute_tau_eps.hpp" +#include "specfem/attenuation/compute_tau_sigma.hpp" +#include "specfem/constants.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/enums.hpp" +#include "specfem/mesh/dim3/materials/materials.hpp" +#include "specfem/setup.hpp" +#include + +namespace specfem::assembly::impl { + +template +struct attenuation_medium : + specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim3> { + + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::attenuation, + specfem::element::dimension_tag::dim3>; + + using view_type = typename base_type::vector_type< + type_real, Kokkos::DefaultExecutionSpace::memory_space>; + + constexpr static int N_SLS = specfem::constants::N_SLS; + + // Host-only per-element scale factors + Kokkos::View h_kappa_scale; + Kokkos::View h_mu_scale; + + // Views: shape [nspec_attn][ngllz][nglly][ngllx][N_SLS] + view_type kappa_relaxation_rate; + view_type::HostMirror h_kappa_relaxation_rate; + view_type mu_relaxation_rate; + view_type::HostMirror h_mu_relaxation_rate; + view_type memory_variable_kappa; + view_type::HostMirror h_memory_variable_kappa; + view_type memory_variable_Rxx; + view_type::HostMirror h_memory_variable_Rxx; + view_type memory_variable_Ryy; + view_type::HostMirror h_memory_variable_Ryy; + view_type memory_variable_Rzz; + view_type::HostMirror h_memory_variable_Rzz; + view_type memory_variable_Rxy; + view_type::HostMirror h_memory_variable_Rxy; + view_type memory_variable_Rxz; + view_type::HostMirror h_memory_variable_Rxz; + view_type memory_variable_Ryz; + view_type::HostMirror h_memory_variable_Ryz; + + // Index mapping: global ispec -> compact attenuation index (-1 if not attenuating) + Kokkos::View + h_attenuation_index_mapping; + Kokkos::View + attenuation_index_mapping; + + attenuation_medium() = default; + + attenuation_medium( + const Kokkos::View elements, + const specfem::assembly::mesh< + specfem::element::dimension_tag::dim3> &mesh, + const specfem::mesh::materials + &materials, + const int ngllz, const int nglly, const int ngllx, + const specfem::units::Hertz fc, const specfem::units::Hertz f0, + const specfem::utilities::Band &band, + const Kokkos::View &tau_sigma) { + + const int nspec_attn = elements.extent(0); + + // 1. Allocate all views + h_kappa_scale = + Kokkos::View( + "kappa_scale", nspec_attn); + h_mu_scale = + Kokkos::View( + "mu_scale", nspec_attn); + kappa_relaxation_rate = + view_type("kappa_relaxation_rate", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_kappa_relaxation_rate = + Kokkos::create_mirror_view(kappa_relaxation_rate); + mu_relaxation_rate = + view_type("mu_relaxation_rate", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_mu_relaxation_rate = + Kokkos::create_mirror_view(mu_relaxation_rate); + memory_variable_kappa = + view_type("mem_kappa", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_kappa = + Kokkos::create_mirror_view(memory_variable_kappa); + memory_variable_Rxx = + view_type("mem_Rxx", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_Rxx = + Kokkos::create_mirror_view(memory_variable_Rxx); + memory_variable_Ryy = + view_type("mem_Ryy", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_Ryy = + Kokkos::create_mirror_view(memory_variable_Ryy); + memory_variable_Rzz = + view_type("mem_Rzz", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_Rzz = + Kokkos::create_mirror_view(memory_variable_Rzz); + memory_variable_Rxy = + view_type("mem_Rxy", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_Rxy = + Kokkos::create_mirror_view(memory_variable_Rxy); + memory_variable_Rxz = + view_type("mem_Rxz", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_Rxz = + Kokkos::create_mirror_view(memory_variable_Rxz); + memory_variable_Ryz = + view_type("mem_Ryz", nspec_attn, ngllz, nglly, ngllx, N_SLS); + h_memory_variable_Ryz = + Kokkos::create_mirror_view(memory_variable_Ryz); + + // Allocate and populate the inverse index mapping (global ispec -> compact index) + h_attenuation_index_mapping = + Kokkos::View( + "h_attenuation_index_mapping", mesh.nspec); + Kokkos::deep_copy(h_attenuation_index_mapping, -1); + for (int i = 0; i < nspec_attn; ++i) { + h_attenuation_index_mapping(elements(i)) = i; + } + attenuation_index_mapping = + Kokkos::View( + "attenuation_index_mapping", mesh.nspec); + Kokkos::deep_copy(attenuation_index_mapping, h_attenuation_index_mapping); + + // Sync zero-initialized device views to host mirrors + copy_to_host(); + + if (nspec_attn == 0) { + return; + } + + // Sanity check input frequencies before looping over elements + if (fc.raw() <= 0 || f0.raw() <= 0) { + throw std::runtime_error( + "Center frequency fc and reference frequency f0 must be positive."); + } + + // 3. Loop over elements + for (int i = 0; i < nspec_attn; ++i) { + const int ispec = elements(i); + + auto material = materials.template get_material< + specfem::element::medium_tag::elastic, PropertyTag, + specfem::element::attenuation_tag::constant_isotropic>(ispec); + + auto computed_values = material.compute_attenuation_properties(fc.raw(), f0.raw(), band, tau_sigma); + auto kappa_props = computed_values.kappa_attenuation_properties; + auto mu_props = computed_values.mu_attenuation_properties; + + // Per-GLL fill + for (int j = 0; j < N_SLS; ++j) { + + // Compute per element relaxation rates + const type_real tauinv_j = 1.0 / tau_sigma(j); + auto kappa_rr_j = kappa_props.beta(j) * tauinv_j / kappa_props.one_minus_sum_beta; + auto mu_rr_j = 2.0 * mu_props.beta(j) * tauinv_j / mu_props.one_minus_sum_beta; + + // Assigning relaxation rates to all GLL points + for (int iz = 0; iz < ngllz; ++iz) { + for (int iy = 0; iy < nglly; ++iy) { + for (int ix = 0; ix < ngllx; ++ix) { + h_kappa_relaxation_rate(i, iz, iy, ix, j) = kappa_rr_j; + h_mu_relaxation_rate(i, iz, iy, ix, j) = mu_rr_j; + } + } + } + } + } + + // 4. Push all host data (kappa/mu_cf filled; memory variables zero) to device + copy_to_device(); + } + + void copy_to_host() { + Kokkos::deep_copy(h_kappa_relaxation_rate, kappa_relaxation_rate); + Kokkos::deep_copy(h_mu_relaxation_rate, mu_relaxation_rate); + Kokkos::deep_copy(h_memory_variable_kappa, memory_variable_kappa); + Kokkos::deep_copy(h_memory_variable_Rxx, memory_variable_Rxx); + Kokkos::deep_copy(h_memory_variable_Ryy, memory_variable_Ryy); + Kokkos::deep_copy(h_memory_variable_Rzz, memory_variable_Rzz); + Kokkos::deep_copy(h_memory_variable_Rxy, memory_variable_Rxy); + Kokkos::deep_copy(h_memory_variable_Rxz, memory_variable_Rxz); + Kokkos::deep_copy(h_memory_variable_Ryz, memory_variable_Ryz); + } + + void copy_to_device() { + Kokkos::deep_copy(kappa_relaxation_rate, h_kappa_relaxation_rate); + Kokkos::deep_copy(mu_relaxation_rate, h_mu_relaxation_rate); + Kokkos::deep_copy(memory_variable_kappa, h_memory_variable_kappa); + Kokkos::deep_copy(memory_variable_Rxx, h_memory_variable_Rxx); + Kokkos::deep_copy(memory_variable_Ryy, h_memory_variable_Ryy); + Kokkos::deep_copy(memory_variable_Rzz, h_memory_variable_Rzz); + Kokkos::deep_copy(memory_variable_Rxy, h_memory_variable_Rxy); + Kokkos::deep_copy(memory_variable_Rxz, h_memory_variable_Rxz); + Kokkos::deep_copy(memory_variable_Ryz, h_memory_variable_Ryz); + } + + /** + * @brief Load attenuation data for a single GLL point from device views + * into a point-local attenuation struct. + * + * Populates relaxation rates and memory variables. The global RK + * coefficients are NOT populated here; they are added by the outer + * load_on_device free function. + */ + template + KOKKOS_INLINE_FUNCTION void load_device_values(const IndexType &index, + PointType &point) const { + const int i = attenuation_index_mapping(index.ispec); + if constexpr (!IndexType::using_simd) { + for (int j = 0; j < N_SLS; ++j) { + point.kappa_relaxation_rate(j) = + kappa_relaxation_rate(i, index.iz, index.iy, index.ix, j); + point.mu_relaxation_rate(j) = + mu_relaxation_rate(i, index.iz, index.iy, index.ix, j); + point.Rxx(j) = memory_variable_Rxx(i, index.iz, index.iy, index.ix, j); + point.Ryy(j) = memory_variable_Ryy(i, index.iz, index.iy, index.ix, j); + point.Rxy(j) = memory_variable_Rxy(i, index.iz, index.iy, index.ix, j); + point.Rxz(j) = memory_variable_Rxz(i, index.iz, index.iy, index.ix, j); + point.Ryz(j) = memory_variable_Ryz(i, index.iz, index.iy, index.ix, j); + point.Rkappa(j) = + memory_variable_kappa(i, index.iz, index.iy, index.ix, j); + } + } else { + using simd = typename PointType::simd; + using mask_type = typename simd::mask_type; + using tag_type = typename simd::tag_type; + const auto mask = index.template get_mask(); + for (int j = 0; j < N_SLS; ++j) { + Kokkos::Experimental::where(mask, point.kappa_relaxation_rate(j)) + .copy_from( + &kappa_relaxation_rate(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.mu_relaxation_rate(j)) + .copy_from( + &mu_relaxation_rate(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxx(j)) + .copy_from( + &memory_variable_Rxx(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Ryy(j)) + .copy_from( + &memory_variable_Ryy(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxy(j)) + .copy_from( + &memory_variable_Rxy(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxz(j)) + .copy_from( + &memory_variable_Rxz(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Ryz(j)) + .copy_from( + &memory_variable_Ryz(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rkappa(j)) + .copy_from( + &memory_variable_kappa(i, index.iz, index.iy, index.ix, j), + tag_type()); + } + } + } + + /** + * @brief Store evolved SLS memory variables from a point-local struct back + * to the device views. + * + * Only the memory variables are written; relaxation rates are + * simulation-lifetime constants and are not written back. + * Note: memory_variable_Rzz = -(Rxx + Ryy) and is not stored in the point + * type; callers must update it separately if needed. + */ + template + KOKKOS_INLINE_FUNCTION void + store_device_values(const IndexType &index, const PointType &point) const { + const int i = attenuation_index_mapping(index.ispec); + if constexpr (!IndexType::using_simd) { + for (int j = 0; j < N_SLS; ++j) { + memory_variable_Rxx(i, index.iz, index.iy, index.ix, j) = point.Rxx(j); + memory_variable_Ryy(i, index.iz, index.iy, index.ix, j) = point.Ryy(j); + memory_variable_Rxy(i, index.iz, index.iy, index.ix, j) = point.Rxy(j); + memory_variable_Rxz(i, index.iz, index.iy, index.ix, j) = point.Rxz(j); + memory_variable_Ryz(i, index.iz, index.iy, index.ix, j) = point.Ryz(j); + memory_variable_kappa(i, index.iz, index.iy, index.ix, j) = + point.Rkappa(j); + } + } else { + using simd = typename PointType::simd; + using mask_type = typename simd::mask_type; + using tag_type = typename simd::tag_type; + const auto mask = index.template get_mask(); + for (int j = 0; j < N_SLS; ++j) { + Kokkos::Experimental::where(mask, point.Rxx(j)) + .copy_to(&memory_variable_Rxx(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Ryy(j)) + .copy_to(&memory_variable_Ryy(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxy(j)) + .copy_to(&memory_variable_Rxy(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rxz(j)) + .copy_to(&memory_variable_Rxz(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Ryz(j)) + .copy_to(&memory_variable_Ryz(i, index.iz, index.iy, index.ix, j), + tag_type()); + Kokkos::Experimental::where(mask, point.Rkappa(j)) + .copy_to( + &memory_variable_kappa(i, index.iz, index.iy, index.ix, j), + tag_type()); + } + } + } +}; + +} // namespace specfem::assembly::impl diff --git a/core/specfem/assembly/attenuation/impl/attenuation_medium/none.tpp b/core/specfem/assembly/attenuation/impl/attenuation_medium/none.tpp new file mode 100644 index 000000000..5a2899443 --- /dev/null +++ b/core/specfem/assembly/attenuation/impl/attenuation_medium/none.tpp @@ -0,0 +1,32 @@ +#pragma once + +#include "specfem/assembly/attenuation/impl/attenuation_medium.hpp" +#include "specfem/enums.hpp" + +namespace specfem::assembly::impl { + +/** + * @brief Empty specialization of attenuation_medium for attenuation_tag::none. + * + * This struct has no data members and all operations are no-ops, ensuring zero + * overhead for element combinations that do not use attenuation. + */ +template +struct attenuation_medium { + attenuation_medium() = default; + void copy_to_host() {} + void copy_to_device() {} + + template + KOKKOS_INLINE_FUNCTION void load_device_values(const IndexType &, + PointType &) const {} + + template + KOKKOS_INLINE_FUNCTION void + store_device_values(const IndexType &, const PointType &) const {} +}; + +} // namespace specfem::assembly::impl diff --git a/core/specfem/assembly/attenuation/load_on_device.hpp b/core/specfem/assembly/attenuation/load_on_device.hpp new file mode 100644 index 000000000..193336c1a --- /dev/null +++ b/core/specfem/assembly/attenuation/load_on_device.hpp @@ -0,0 +1,69 @@ +#pragma once + +#include "specfem/assembly/attenuation.hpp" +#include "specfem/constants.hpp" +#include "specfem/enums.hpp" +#include "specfem/setup.hpp" +#include + +namespace specfem::assembly { + +/** + * @defgroup AttenuationDataAccess Attenuation Data Access Functions + */ + +/** + * @brief Load attenuation state for a GLL point from the assembly attenuation + * container into a point-local attenuation struct (device kernel). + * + * Populates all fields of @p point: + * - kappa and mu relaxation rates (per-element constants read from device + * views inside the matching @c attenuation_medium), + * - the current SLS memory variables (e.g. Rxx, Rxz, Rkappa for dim2), and + * - the global Runge-Kutta coefficients (alpha/beta/gamma, broadcast from the + * @c Attenuation device views to all active SIMD lanes). + * + * @tparam PointAttenuationType Point attenuation struct + * (specfem::point::attenuation<...>). + * Must have @c medium_tag, @c property_tag, and + * @c attenuation_tag == constant_isotropic. + * @tparam IndexType Point index struct + * (specfem::point::index<...>). + * @c IndexType::using_simd must match + * @c PointAttenuationType::using_simd. + * + * @param lcoord GLL point index (global ispec, iz, [iy,] ix). + * @param attenuation Assembly attenuation container (device-accessible). + * @param point Output: populated point attenuation struct. + * + * @ingroup AttenuationDataAccess + */ +template = 0> +KOKKOS_FORCEINLINE_FUNCTION void load_on_device( + const IndexType &lcoord, + const specfem::assembly::Attenuation + &attenuation, + PointAttenuationType &point) { + + constexpr auto MediumTag = PointAttenuationType::medium_tag; + constexpr auto PropertyTag = PointAttenuationType::property_tag; + constexpr int N_SLS = PointAttenuationType::N_SLS; + + // Load per-element data: relaxation rates and memory variables + attenuation.template get_container() + .load_device_values(lcoord, point); + + // Broadcast global RK coefficients (identical for every element in a run) + for (int j = 0; j < N_SLS; ++j) { + point.alpha_rk(j) = attenuation.d_alpha_rk(j); + point.beta_rk(j) = attenuation.d_beta_rk(j); + point.gamma_rk(j) = attenuation.d_gamma_rk(j); + } +} + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/attenuation/store_on_device.hpp b/core/specfem/assembly/attenuation/store_on_device.hpp new file mode 100644 index 000000000..408ae1640 --- /dev/null +++ b/core/specfem/assembly/attenuation/store_on_device.hpp @@ -0,0 +1,56 @@ +#pragma once + +#include "specfem/assembly/attenuation.hpp" +#include "specfem/enums.hpp" +#include + +namespace specfem::assembly { + +/** + * @brief Store evolved SLS memory variables from a point-local attenuation + * struct back to the assembly attenuation container (device kernel). + * + * Only the memory variables (e.g. Rxx, Rxz, Rkappa for dim2; additionally + * Ryy, Rxy, Rxz, Ryz for dim3) are written back. The relaxation rates and + * Runge-Kutta coefficients are simulation-lifetime constants and are never + * written back by this function. + * + * @note For dim3, @c memory_variable_Rzz = -(Rxx + Ryy) is not stored in the + * point type and must be updated separately by the caller if required. + * + * @tparam PointAttenuationType Point attenuation struct + * (specfem::point::attenuation<...>). + * Must have @c medium_tag, @c property_tag, and + * @c attenuation_tag == constant_isotropic. + * @tparam IndexType Point index struct + * (specfem::point::index<...>). + * @c IndexType::using_simd must match + * @c PointAttenuationType::using_simd. + * + * @param index GLL point index (global ispec, iz, [iy,] ix). + * @param attenuation Assembly attenuation container (device-accessible). + * @param point Input: point attenuation struct with updated memory vars. + * + * @ingroup AttenuationDataAccess + */ +template = 0> +KOKKOS_FORCEINLINE_FUNCTION void store_on_device( + const IndexType &index, + const specfem::assembly::Attenuation + &attenuation, + const PointAttenuationType &point) { + + constexpr auto MediumTag = PointAttenuationType::medium_tag; + constexpr auto PropertyTag = PointAttenuationType::property_tag; + + // Store memory variables back to device views + attenuation.template get_container() + .store_device_values(index, point); +} + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/boundaries/boundaries.hpp b/core/specfem/assembly/boundaries/boundaries.hpp new file mode 100644 index 000000000..7aaf18591 --- /dev/null +++ b/core/specfem/assembly/boundaries/boundaries.hpp @@ -0,0 +1,32 @@ +#pragma once + +#include "specfem/enums.hpp" + +namespace specfem::assembly::boundaries_impl { + +/** + * @brief Primary template declaration for acoustic free surface boundary + * conditions. + */ +template +struct acoustic_free_surface; + +/** + * @brief Primary template declaration for Stacey boundary conditions. + */ +template struct stacey; + +} // namespace specfem::assembly::boundaries_impl + +namespace specfem::assembly { + +/** + * @brief Data container for boundary condition information at every quadrature + * point on the boundary. + * + * Primary template declaration. Dimension-specific partial specializations are + * defined in boundaries/dim2/boundaries.hpp and boundaries/dim3/boundaries.hpp. + */ +template class boundaries; + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/element_types/dim3/element_types.cpp b/core/specfem/assembly/element_types/dim3/element_types.cpp index 3a057ea9d..4483e70f1 100644 --- a/core/specfem/assembly/element_types/dim3/element_types.cpp +++ b/core/specfem/assembly/element_types/dim3/element_types.cpp @@ -44,7 +44,7 @@ specfem::assembly::element_types:: FOR_EACH_IN_PRODUCT( (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC, ACOUSTIC), - PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE)), + PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)), CAPTURE(material_elements, h_material_elements) { int count = 0; int index = 0; @@ -138,15 +138,15 @@ Kokkos::View specfem::assembly:: const specfem::element::property_tag property_tag, const specfem::element::attenuation_tag attenuation_tag) const { - FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC, ACOUSTIC), - PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE)), - CAPTURE(h_material_elements) { - if (_medium_tag_ == medium_tag && - _property_tag_ == property_tag && - _attenuation_tag_ == attenuation_tag) { - return _h_material_elements_; - } - }) + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC, ACOUSTIC), + PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)), + CAPTURE(h_material_elements) { + if (_medium_tag_ == medium_tag && _property_tag_ == property_tag && + _attenuation_tag_ == attenuation_tag) { + return _h_material_elements_; + } + }) throw std::runtime_error("Medium tag or property tag not found"); } @@ -158,15 +158,15 @@ specfem::assembly::element_types:: const specfem::element::property_tag property_tag, const specfem::element::attenuation_tag attenuation_tag) const { - FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC, ACOUSTIC), - PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE)), - CAPTURE(material_elements) { - if (_medium_tag_ == medium_tag && - _property_tag_ == property_tag && - _attenuation_tag_ == attenuation_tag) { - return _material_elements_; - } - }) + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC, ACOUSTIC), + PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)), + CAPTURE(material_elements) { + if (_medium_tag_ == medium_tag && _property_tag_ == property_tag && + _attenuation_tag_ == attenuation_tag) { + return _material_elements_; + } + }) throw std::runtime_error("Medium tag or property tag not found"); } diff --git a/core/specfem/assembly/field_derivative_storage.hpp b/core/specfem/assembly/field_derivative_storage.hpp new file mode 100644 index 000000000..0b27d0f7b --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage.hpp @@ -0,0 +1,13 @@ +#pragma once + +#include "specfem/enums.hpp" + +namespace specfem::assembly { + +template +struct FieldDerivativeStorage; + +} // namespace specfem::assembly + +#include "specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.hpp" +#include "specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.hpp" diff --git a/core/specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.cpp b/core/specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.cpp new file mode 100644 index 000000000..7f8ece4bd --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.cpp @@ -0,0 +1,46 @@ +#include "specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.hpp" +#include "specfem/assembly/element_types.hpp" +#include "specfem/macros.hpp" +#include + +specfem::assembly::FieldDerivativeStorage< + specfem::element::dimension_tag::dim2>:: + FieldDerivativeStorage( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim2> &element_types, + const int nspec_global, const int ngllz, const int ngllx) { + + // Initialize storage for all CONSTANT_ISOTROPIC combinations. + // NONE combinations are default-constructed (empty, no-op). + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV)), CAPTURE(fd_medium) { + // Select elements in this medium that require field derivative storage + // (e.g. are attenuating). + std::vector element_indices; + + auto elements = element_types.get_elements_on_host(_medium_tag_); + + for (int element_index = 0; element_index < elements.extent(0); + ++element_index) { + + const int ispec = elements(element_index); + + // This expresion could be + bool requires_storage = (element_types.attenuation_tags(ispec) != + specfem::element::attenuation_tag::none); + + // Store the global + if (requires_storage) { + element_indices.push_back(ispec); + } + } + + const size_t count = element_indices.size(); + Kokkos::View subset_elements( + "subset_elements", count); + for (size_t i = 0; i < count; ++i) + subset_elements(i) = element_indices[i]; + + _fd_medium_ = { subset_elements, nspec_global, ngllz, ngllx }; + }) +} diff --git a/core/specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.hpp b/core/specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.hpp new file mode 100644 index 000000000..aedb309b5 --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/dim2/field_derivative_storage.hpp @@ -0,0 +1,97 @@ +#pragma once + +#include "specfem/assembly/element_types.hpp" +#include "specfem/assembly/field_derivative_storage.hpp" +#include "specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros.hpp" +#include "specfem/setup.hpp" +#include +#include + +namespace specfem::assembly { + +/** + * @brief 2D per-GLL-point field derivative storage for attenuation strain + * bookkeeping. + * + * Holds one field_derivative_medium per medium tag. But selects whether to + * store field_derivative values based on the, e.g., attenuation tag: for + * attenuation_none, the + * + * Follows the same pattern as specfem::assembly::properties. + */ +template <> +struct FieldDerivativeStorage + : specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::field_derivatives, + specfem::element::dimension_tag::dim2> { + + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::field_derivatives, + specfem::element::dimension_tag::dim2>; + + constexpr static auto dimension_tag = specfem::element::dimension_tag::dim2; + + // One field_derivative_medium per medium, but only some have storage + // depending on the attenuation tag. For combination. Attenuation_none → empty + // struct. Constant_isotropic → compact storage. + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV)), + DECLARE(((specfem::assembly::field_derivative_storage:: + impl::field_derivative_medium, + (_DIMENSION_TAG_, _MEDIUM_TAG_)), + fd_medium))) + + FieldDerivativeStorage() = default; + + /** + * @brief Construct and allocate storage for all non-none element + * combinations. + * + * @param element_types Element type information (provides per-combination + * element index lists). + * @param nspec_global Total spectral elements (for ispec_to_compact + * mapping). + * @param ngllz GLL points in z-direction. + * @param ngllx GLL points in x-direction. + */ + FieldDerivativeStorage( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim2> &element_types, + const int nspec_global, const int ngllz, const int ngllx); + + /** + * @brief Access the field_derivative_medium for a given tag combination. + */ + template + KOKKOS_INLINE_FUNCTION constexpr specfem::assembly::field_derivative_storage:: + impl::field_derivative_medium const & + get_container() const { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV)), + CAPTURE(fd_medium) { + if constexpr (_medium_tag_ == MediumTag) { + return _fd_medium_; + } + }) + Kokkos::abort( + "Invalid tag combination in FieldDerivativeStorage::get_medium"); + SUPPRESS_UNREACHABLE(return {};) + } + + void copy_to_host() { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV)), + CAPTURE(fd_medium) { _fd_medium_.copy_to_host(); }) + } + + void copy_to_device() { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM2), MEDIUM_TAG(ELASTIC_PSV)), + CAPTURE(fd_medium) { _fd_medium_.copy_to_device(); }) + } +}; + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.cpp b/core/specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.cpp new file mode 100644 index 000000000..17b1cf317 --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.cpp @@ -0,0 +1,45 @@ +#include "specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.hpp" +#include "specfem/assembly/element_types.hpp" +#include "specfem/macros.hpp" +#include + +specfem::assembly::FieldDerivativeStorage< + specfem::element::dimension_tag::dim3>:: + FieldDerivativeStorage( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim3> &element_types, + const int nspec_global, const int ngllz, const int nglly, + const int ngllx) { + + FOR_EACH_IN_PRODUCT( + (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC)), CAPTURE(fd_medium) { + // Select elements in this medium that require field derivative storage + // (i.e. are attenuating). + std::vector element_indices; + + auto elements = element_types.get_elements_on_host(_medium_tag_); + + for (int element_index = 0; element_index < elements.extent(0); + ++element_index) { + + const int ispec = elements(element_index); + + // This expresion could be + bool requires_storage = (element_types.attenuation_tags(ispec) != + specfem::element::attenuation_tag::none); + + // Store the global + if (requires_storage) { + element_indices.push_back(ispec); + } + } + + const size_t count = element_indices.size(); + Kokkos::View subset_elements( + "subset_elements", count); + for (size_t i = 0; i < count; ++i) + subset_elements(i) = element_indices[i]; + + _fd_medium_ = { subset_elements, nspec_global, ngllz, nglly, ngllx }; + }) +} diff --git a/core/specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.hpp b/core/specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.hpp new file mode 100644 index 000000000..9112137ed --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/dim3/field_derivative_storage.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include "specfem/assembly/element_types.hpp" +#include "specfem/assembly/field_derivative_storage.hpp" +#include "specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp" +#include "specfem/assembly/mesh.hpp" +#include "specfem/data_access/container.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros.hpp" +#include "specfem/setup.hpp" +#include +#include + +namespace specfem::assembly { + +/** + * @brief 3D per-GLL-point field derivative storage for attenuation field + * derivatives bookkeeping. + * + * Follows the same pattern as the 2D specialization. For attenuation_none + * combinations the sub-struct is empty (zero overhead). + */ +template <> +struct FieldDerivativeStorage + : specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::field_derivatives, + specfem::element::dimension_tag::dim3> { + + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::field_derivatives, + specfem::element::dimension_tag::dim3>; + + constexpr static auto dimension_tag = specfem::element::dimension_tag::dim3; + + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC)), + DECLARE(((specfem::assembly::field_derivative_storage:: + impl::field_derivative_medium, + (_DIMENSION_TAG_, _MEDIUM_TAG_)), + fd_medium))) + + FieldDerivativeStorage() = default; + + FieldDerivativeStorage( + const specfem::assembly::element_types< + specfem::element::dimension_tag::dim3> &element_types, + const int nspec_global, const int ngllz, const int nglly, + const int ngllx); + + template + KOKKOS_INLINE_FUNCTION constexpr specfem::assembly::field_derivative_storage:: + impl::field_derivative_medium const & + get_container() const { + + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC)), + CAPTURE(fd_medium) { + if constexpr (_medium_tag_ == MediumTag) { + return _fd_medium_; + } + }) + Kokkos::abort( + "Invalid tag combination in FieldDerivativeStorage::get_container"); + + SUPPRESS_UNREACHABLE(return {};) + } + + void copy_to_host() { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC)), + CAPTURE(fd_medium) { _fd_medium_.copy_to_host(); }) + } + + void copy_to_device() { + FOR_EACH_IN_PRODUCT((DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC)), + CAPTURE(fd_medium) { _fd_medium_.copy_to_device(); }) + } +}; + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp b/core/specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp new file mode 100644 index 000000000..4c5de63ab --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp @@ -0,0 +1,272 @@ +#pragma once + +#include "specfem/data_access/container.hpp" +#include "specfem/element.hpp" +#include "specfem/enums.hpp" +#include "specfem/setup.hpp" +#include +#include + +namespace specfem::assembly::field_derivative_storage::impl { + +/** + * @brief Primary template for field derivative storage at a set of GLL points. + * + * Stores the du (strain) tensor per GLL point for elements belonging to a given + * DimensionTag, MediumTag combination. Used to record + * the field derivatives from the previous time step for attenuation strain + * calculations. + * + * Storage is compact: only elements with this specific tag combination are + * stored, indexed via a global ispec → compact index mapping. + * + * This primary template covers all AttenuationTag values except + * attenuation_none. For attenuation_none, an empty specialization below + * provides zero overhead. + * + * @tparam DimensionTag Dimension of the elements + * @tparam MediumTag Medium of the elements + */ +template +struct field_derivative_medium + : specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::field_derivatives, + DimensionTag> { + + using base_type = specfem::data_access::Container< + specfem::data_access::ContainerType::domain, + specfem::data_access::DataClassType::field_derivatives, DimensionTag>; + + constexpr static auto dimension_tag = DimensionTag; + constexpr static auto medium_tag = MediumTag; + + /** + * @name Compile-time constants + */ + ///@{ + static constexpr int components = + specfem::element::attributes::components; + static constexpr int num_dimensions = + specfem::element::attributes::dimension; + ///@} + + using view_type = typename base_type::template tensor_type< + type_real, Kokkos::DefaultExecutionSpace::memory_space>; + + /// Integer index view for the compact → global ispec mapping + using index_view_type = Kokkos::View; + + /// Compact storage: shape + /// [nspec_attn][ngllz][...][components][num_dimensions] + view_type du; + typename view_type::HostMirror h_du; + + /// Maps global ispec → compact index (or -1 for non-matching elements) + index_view_type ispec_to_compact; + typename index_view_type::HostMirror h_ispec_to_compact; + + field_derivative_medium() = default; + + /** + * @brief Construct storage for the given element list. + * + * @param elements Host view of global ispec indices belonging to this + * medium + * @param nspec_global Total number of spectral elements in the mesh (for + * the inverse mapping view). + * @param ngllz Number of GLL points in z-direction. + * @param ngllx Number of GLL points in x-direction. + */ + template < + specfem::element::dimension_tag D = DimensionTag, + std::enable_if_t = 0> + field_derivative_medium( + const Kokkos::View &elements, + const int nspec_global, const int ngllz, const int ngllx) { + + const int nspec_attn = static_cast(elements.extent(0)); + + du = view_type("field_derivative_storage", nspec_attn, ngllz, ngllx, + components, num_dimensions); + h_du = typename view_type::HostMirror("h_field_derivative_storage", + nspec_attn, ngllz, ngllx, components, + num_dimensions); + Kokkos::deep_copy(du, static_cast(0)); + + // Build ispec → compact-index inverse mapping (initialized to -1) + ispec_to_compact = + index_view_type("ispec_to_compact_field_deriv", nspec_global); + h_ispec_to_compact = Kokkos::create_mirror_view(ispec_to_compact); + Kokkos::deep_copy(h_ispec_to_compact, -1); + + for (int i = 0; i < nspec_attn; ++i) { + h_ispec_to_compact(elements(i)) = i; + } + + copy_to_host(); // sync h_du from zeroed device view + copy_to_device(); // push inverse mapping to device + } + + /** + * @brief Construct storage for the given element list. + * + * @param elements Host view of global ispec indices belonging to this + * medium. + * @param nspec_global Total number of spectral elements in the mesh (for + * the inverse mapping view). + * @param ngllz Number of GLL points in z-direction. + * @param nglly Number of GLL points in y-direction + * @param ngllx Number of GLL points in x-direction. + */ + template < + specfem::element::dimension_tag D = DimensionTag, + std::enable_if_t = 0> + field_derivative_medium( + const Kokkos::View &elements, + const int nspec_global, const int ngllz, const int nglly, + const int ngllx) { + + const int nspec_attn = static_cast(elements.extent(0)); + + du = view_type("field_derivative_storage", nspec_attn, ngllz, nglly, ngllx, + components, num_dimensions); + h_du = typename view_type::HostMirror("h_field_derivative_storage", + nspec_attn, ngllz, nglly, ngllx, + components, num_dimensions); + + Kokkos::deep_copy(du, static_cast(0)); + + // Build ispec → compact-index inverse mapping (initialized to -1) + ispec_to_compact = + index_view_type("ispec_to_compact_field_deriv", nspec_global); + h_ispec_to_compact = Kokkos::create_mirror_view(ispec_to_compact); + Kokkos::deep_copy(h_ispec_to_compact, -1); + for (int i = 0; i < nspec_attn; ++i) { + h_ispec_to_compact(elements(i)) = i; + } + + copy_to_host(); // sync h_du from zeroed device view + copy_to_device(); // push inverse mapping to device + } + + void copy_to_host() { + Kokkos::deep_copy(h_du, du); + Kokkos::deep_copy(h_ispec_to_compact, ispec_to_compact); + } + + void copy_to_device() { + Kokkos::deep_copy(du, h_du); + Kokkos::deep_copy(ispec_to_compact, h_ispec_to_compact); + } + + /** + * @brief Load field derivatives for a non-SIMD GLL point from compact + * storage. + */ + template = 0> + KOKKOS_FORCEINLINE_FUNCTION void + load_device_values(const IndexType &index, PointFDType &point_fd) const { + const int i = ispec_to_compact(index.ispec); + for (int ic = 0; ic < PointFDType::components; ++ic) { + for (int id = 0; id < PointFDType::num_dimensions; ++id) { + if constexpr (DimensionTag == specfem::element::dimension_tag::dim2) { + point_fd.du[ic][id] = du(i, index.iz, index.ix, ic, id); + } else { + point_fd.du[ic][id] = du(i, index.iz, index.iy, index.ix, ic, id); + } + } + } + } + + /** + * @brief Load field derivatives for a SIMD GLL point from compact storage. + * + * Uses copy_from to load SIMD_WIDTH consecutive compact-index values + * starting at the compact index of the first lane, matching the + * jacobian_matrix SIMD pattern. + */ + template = 0> + KOKKOS_FORCEINLINE_FUNCTION void + load_device_values(const IndexType &index, PointFDType &point_fd) const { + using simd_type = typename PointFDType::simd; + using mask_type = typename simd_type::mask_type; + using tag_type = typename simd_type::tag_type; + + const int i = ispec_to_compact(index.ispec); + const auto mask = index.template get_mask(); + + for (int ic = 0; ic < PointFDType::components; ++ic) { + for (int id = 0; id < PointFDType::num_dimensions; ++id) { + if constexpr (DimensionTag == specfem::element::dimension_tag::dim2) { + Kokkos::Experimental::where(mask, point_fd.du[ic][id]) + .copy_from(&du(i, index.iz, index.ix, ic, id), tag_type()); + } else { + Kokkos::Experimental::where(mask, point_fd.du[ic][id]) + .copy_from(&du(i, index.iz, index.iy, index.ix, ic, id), + tag_type()); + } + } + } + } + + /** + * @brief Store field derivatives for a non-SIMD GLL point into compact + * storage. + */ + template = 0> + KOKKOS_FORCEINLINE_FUNCTION void + store_device_values(const IndexType &index, const PointFDType &point_fd) { + const int i = ispec_to_compact(index.ispec); + for (int ic = 0; ic < PointFDType::components; ++ic) { + for (int id = 0; id < PointFDType::num_dimensions; ++id) { + if constexpr (DimensionTag == specfem::element::dimension_tag::dim2) { + du(i, index.iz, index.ix, ic, id) = point_fd.du[ic][id]; + } else { + du(i, index.iz, index.iy, index.ix, ic, id) = point_fd.du[ic][id]; + } + } + } + } + + /** + * @brief Store field derivatives for a SIMD GLL point into compact storage. + * + * Uses copy_to to scatter SIMD_WIDTH consecutive compact-index values + * starting at the compact index of the first lane, matching the + * jacobian_matrix SIMD pattern. + */ + template = 0> + KOKKOS_FORCEINLINE_FUNCTION void + store_device_values(const IndexType &index, const PointFDType &point_fd) { + using simd_type = typename PointFDType::simd; + using mask_type = typename simd_type::mask_type; + using tag_type = typename simd_type::tag_type; + + const int i = ispec_to_compact(index.ispec); + const auto mask = index.template get_mask(); + + for (int ic = 0; ic < PointFDType::components; ++ic) { + for (int id = 0; id < PointFDType::num_dimensions; ++id) { + if constexpr (DimensionTag == specfem::element::dimension_tag::dim2) { + Kokkos::Experimental::where(mask, point_fd.du[ic][id]) + .copy_to(&du(i, index.iz, index.ix, ic, id), tag_type()); + } else { + Kokkos::Experimental::where(mask, point_fd.du[ic][id]) + .copy_to(&du(i, index.iz, index.iy, index.ix, ic, id), + tag_type()); + } + } + } + } +}; + +} // namespace specfem::assembly::field_derivative_storage::impl diff --git a/core/specfem/assembly/field_derivative_storage/load_on_device.hpp b/core/specfem/assembly/field_derivative_storage/load_on_device.hpp new file mode 100644 index 000000000..78c38ebea --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/load_on_device.hpp @@ -0,0 +1,80 @@ +#pragma once + +#include "specfem/assembly/field_derivative_storage.hpp" +#include "specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp" +#include "specfem/data_access/check_compatibility.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros.hpp" +#include "specfem/setup.hpp" +#include + +namespace specfem::assembly { + +/** + * @defgroup FieldDerivativeStorageDataAccess Field Derivative Storage Data + * Access Functions + * + */ + +/** + * @brief Load field derivatives for a GLL point from the assembly + * FieldDerivativeStorage container (device kernel). + * + * Dispatches to the matching field_derivative_medium via get_container() + * and delegates to its load_device_values() method. + * + * If @p point_fd is specfem::point::null_field_derivatives this function + * returns immediately — no storage exists for this element combination. + * + * @ingroup FieldDerivativeStorageDataAccess + * + * @tparam IndexType Index type + * (specfem::point::index<...>) + * @tparam FieldDerivativesContainerType Assembly FieldDerivativeStorage + * @tparam PointFieldDerivativesType specfem::point::field_derivatives<...> + * or + * specfem::point::null_field_derivatives + * + * @param index GLL point index (global ispec, iz, [iy,] ix). + * @param container Assembly field-derivative storage container + * (device-accessible). + * @param point_fd Output: populated point field-derivatives struct, or + * null_field_derivatives for a no-op. + */ +template ::value && + specfem::data_access::is_field_derivatives< + FieldDerivativesContainerType>::value && + specfem::data_access::is_field_derivatives< + PointFieldDerivativesType>::value> int> += 0 > + KOKKOS_FORCEINLINE_FUNCTION void load_on_device( + const IndexType &index, const FieldDerivativesContainerType &container, + const impl::field_derivatives_type &point_fd) { + + auto compatibility = specfem::data_access::check_compatibility< + IndexType, FieldDerivativesContainerType, PointFieldDerivativesType>(); + + // Incompatible types; this likely indicates a programming error in the + // caller. Trigger a compile-time error with a helpful message. + static_assert(compatibility::value, + "Incompatible types for field derivative loading"); + + constexpr auto medium_tag = + PointFieldDerivativesType::field_derivatives_medium_type::medium_tag; + + container.template get_container().load_device_values(index, + point_fd); +} + +template +KOKKOS_FORCEINLINE_FUNCTION void +load_on_device(const IndexType &index, + const FieldDerivativesContainerType &container, + specfem::point::null_field_derivatives_type point_fd) { + // No-op for null_field_derivatives (no storage allocated) +} + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/field_derivative_storage/store_on_device.hpp b/core/specfem/assembly/field_derivative_storage/store_on_device.hpp new file mode 100644 index 000000000..facc9972e --- /dev/null +++ b/core/specfem/assembly/field_derivative_storage/store_on_device.hpp @@ -0,0 +1,74 @@ +#pragma once + +#include "specfem/assembly/field_derivative_storage.hpp" +#include "specfem/assembly/field_derivative_storage/impl/field_derivative_medium.hpp" +#include "specfem/data_access/check_compatibility.hpp" +#include "specfem/enums.hpp" +#include "specfem/macros.hpp" +#include "specfem/setup.hpp" +#include + +namespace specfem::assembly { + +/** + * @brief Store field derivatives for a GLL point into the assembly + * FieldDerivativeStorage container (device kernel). + * + * Dispatches to the matching field_derivative_medium via get_container() + * and delegates to its store_device_values() method. + * + * If @p point_fd is specfem::point::null_field_derivatives this function + * returns immediately — no storage exists for this element combination. + * + * @ingroup FieldDerivativeStorageDataAccess + * + * @tparam IndexType Index type + * (specfem::point::index<...>) + * @tparam FieldDerivativesContainerType Assembly FieldDerivativeStorage + * @tparam PointFieldDerivativesType specfem::point::field_derivatives<...> + * or + * specfem::point::null_field_derivatives + * + * @param index GLL point index (global ispec, iz, [iy,] ix). + * @param container Assembly field-derivative storage container + * (device-accessible). + * @param point_fd Input: point field-derivatives struct with values to store, + * or null_field_derivatives for a no-op. + */ +template ::value && + specfem::data_access::is_field_derivatives< + FieldDerivativesContainerType>::value && + specfem::data_access::is_field_derivatives< + PointFieldDerivativesType>::value, + int> = 0> +KOKKOS_FORCEINLINE_FUNCTION void +store_on_device(const IndexType &index, + const FieldDerivativesContainerType &container, + const PointFieldDerivativesType &point_fd) { + + auto compatibility = specfem::data_access::check_compatibility< + IndexType, FieldDerivativesContainerType, PointFieldDerivativesType>(); + // Incompatible types; this likely indicates a programming error in the + // caller. Trigger a compile-time error with a helpful message. + static_assert(compatibility::value, + "Incompatible types for field derivative storing."); + + constexpr auto medium_tag = + PointFieldDerivativesType::field_derivatives_medium_type::medium_tag; + + container.template get_container().store_device_values(index, + point_fd); +} + +template +KOKKOS_FORCEINLINE_FUNCTION void +store_on_device(const IndexType &index, + const FieldDerivativesContainerType &container, + specfem::point::null_field_derivatives_type point_fd) { + // No-op for null_field_derivatives (no storage allocated) +} + +} // namespace specfem::assembly diff --git a/core/specfem/assembly/impl/domain_properties.hpp b/core/specfem/assembly/impl/domain_properties.hpp index fc89498cd..359cd19f8 100644 --- a/core/specfem/assembly/impl/domain_properties.hpp +++ b/core/specfem/assembly/impl/domain_properties.hpp @@ -94,6 +94,7 @@ struct domain_properties domain_properties( const Kokkos::View elements, const specfem::assembly::mesh &mesh, const int ngllz, @@ -101,7 +102,9 @@ struct domain_properties - property_index_mapping); + property_index_mapping, + std::integral_constant); /// Device value access disabled for this container type template @@ -164,13 +167,16 @@ struct domain_properties domain_properties( const Kokkos::View elements, const int nspec, const int ngllz, const int nglly, const int ngllx, const specfem::mesh::materials &materials, const Kokkos::View - property_index_mapping); + property_index_mapping, + std::integral_constant); /// Device value access disabled for this container type template diff --git a/core/specfem/assembly/impl/domain_properties.tpp b/core/specfem/assembly/impl/domain_properties.tpp index af8c4a664..fa8028034 100644 --- a/core/specfem/assembly/impl/domain_properties.tpp +++ b/core/specfem/assembly/impl/domain_properties.tpp @@ -4,6 +4,7 @@ template +template specfem::assembly::impl::domain_properties:: domain_properties( @@ -14,7 +15,8 @@ specfem::assembly::impl::domain_properties - property_index_mapping) + property_index_mapping, + std::integral_constant) : base_type(elements.extent(0), ngllz, ngllx) { const int nelement = elements.extent(0); @@ -29,7 +31,7 @@ specfem::assembly::impl::domain_properties(mesh_ispec); + AttenuationTag>(mesh_ispec); // Assign the material property to the property container auto point_property = material.get_properties(); @@ -51,6 +53,7 @@ specfem::assembly::impl::domain_properties +template specfem::assembly::impl::domain_properties:: domain_properties( @@ -59,7 +62,8 @@ specfem::assembly::impl::domain_properties &materials, const Kokkos::View - property_index_mapping) + property_index_mapping, + std::integral_constant) : base_type(elements.extent(0), ngllz, nglly, ngllx) { const int nelement = elements.extent(0); @@ -77,7 +81,8 @@ specfem::assembly::impl::domain_properties(ispec); + AttenuationTag>(ispec); + this->store_host_values( specfem::point::index(count, iz, iy, ix), point_property); diff --git a/core/specfem/assembly/impl/value_containers/dim2/value_containers.hpp b/core/specfem/assembly/impl/value_containers/dim2/value_containers.hpp index 08d7404be..29c0c6215 100644 --- a/core/specfem/assembly/impl/value_containers/dim2/value_containers.hpp +++ b/core/specfem/assembly/impl/value_containers/dim2/value_containers.hpp @@ -42,7 +42,7 @@ struct value_containers:: MEDIUM_TAG(ELASTIC_PSV, ELASTIC_SH, ACOUSTIC, POROELASTIC, ELASTIC_PSV_T), PROPERTY_TAG(ISOTROPIC, ANISOTROPIC, ISOTROPIC_COSSERAT), - ATTENUATION_TAG(NONE)), + ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)), CAPTURE(value) { _value_ = specfem::assembly::impl::domain_properties< _dimension_tag_, _medium_tag_, _property_tag_>( element_types.get_elements_on_host(_medium_tag_, _property_tag_, _attenuation_tag_), mesh, ngllz, ngllx, materials, has_gll_model, - h_property_index_mapping); + h_property_index_mapping, + std::integral_constant{}); }) Kokkos::deep_copy(property_index_mapping, h_property_index_mapping); diff --git a/core/specfem/assembly/properties/dim3/properties.cpp b/core/specfem/assembly/properties/dim3/properties.cpp index bb7d58833..3aed1d0d9 100644 --- a/core/specfem/assembly/properties/dim3/properties.cpp +++ b/core/specfem/assembly/properties/dim3/properties.cpp @@ -26,13 +26,15 @@ specfem::assembly::properties:: FOR_EACH_IN_PRODUCT( (DIMENSION_TAG(DIM3), MEDIUM_TAG(ELASTIC, ACOUSTIC), - PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE)), + PROPERTY_TAG(ISOTROPIC), ATTENUATION_TAG(NONE, CONSTANT_ISOTROPIC)), CAPTURE(value) { _value_ = specfem::assembly::impl::domain_properties< _dimension_tag_, _medium_tag_, _property_tag_>( element_types.get_elements_on_host(_medium_tag_, _property_tag_, _attenuation_tag_), - nspec, ngllz, nglly, ngllx, materials, h_property_index_mapping); + nspec, ngllz, nglly, ngllx, materials, h_property_index_mapping, + std::integral_constant{}); }) Kokkos::deep_copy(property_index_mapping, h_property_index_mapping); diff --git a/core/specfem/attenuation/attenuation.cpp b/core/specfem/attenuation/attenuation.cpp index 6eb643adf..eedc870b6 100644 --- a/core/specfem/attenuation/attenuation.cpp +++ b/core/specfem/attenuation/attenuation.cpp @@ -5,19 +5,6 @@ #include "specfem/attenuation/compute_tau_eps.hpp" #include "specfem/attenuation/compute_tau_sigma.hpp" -template Kokkos::View -specfem::attenuation::compute_tau_sigma( - const type_real, const type_real); - -template Kokkos::View - specfem::attenuation::compute_tau_eps( - type_real, - Kokkos::View, - type_real, type_real); - template specfem::attenuation::AttenuationPropertyValues< specfem::constants::N_SLS> specfem::attenuation::get_attenuation_property_values< diff --git a/core/specfem/attenuation/compute_band.hpp b/core/specfem/attenuation/compute_band.hpp index bfb305659..059d0b60b 100644 --- a/core/specfem/attenuation/compute_band.hpp +++ b/core/specfem/attenuation/compute_band.hpp @@ -32,10 +32,10 @@ namespace attenuation { * * @tparam N_SLS Number of standard linear solids (must be 2–5) * @param min_resolved_period Minimum period resolved by the mesh (s) - * @return Band with @c min and @c max + * @return Band in Hertz with @c min and @c max */ template -specfem::utilities::Band +specfem::utilities::Band compute_band(specfem::units::Seconds min_resolved_period); } // namespace attenuation diff --git a/core/specfem/attenuation/compute_band.tpp b/core/specfem/attenuation/compute_band.tpp index 95f59aa18..e8172453e 100644 --- a/core/specfem/attenuation/compute_band.tpp +++ b/core/specfem/attenuation/compute_band.tpp @@ -8,7 +8,7 @@ namespace specfem { namespace attenuation { template -specfem::utilities::Band +specfem::utilities::Band compute_band(const specfem::units::Seconds min_resolved_period) { static_assert(N_SLS >= 2 && N_SLS <= 5, "N_SLS must be between 2 and 5 (inclusive)"); @@ -20,9 +20,7 @@ compute_band(const specfem::units::Seconds min_resolved_period) { const specfem::units::Seconds max_period = min_resolved_period * static_cast(std::pow(10.0, theta[N_SLS])); - return specfem::utilities::Band{ - min_resolved_period, max_period - }; + return specfem::utilities::Band(1.0/max_period, 1.0/min_resolved_period); } } // namespace attenuation diff --git a/core/specfem/attenuation/compute_tau_eps.hpp b/core/specfem/attenuation/compute_tau_eps.hpp index ad5f52ddf..e1b083ceb 100644 --- a/core/specfem/attenuation/compute_tau_eps.hpp +++ b/core/specfem/attenuation/compute_tau_eps.hpp @@ -1,10 +1,10 @@ #pragma once -#include "compute_tau_sigma.hpp" #include "maxwell.hpp" #include "specfem/constants.hpp" #include "specfem/optimization.hpp" #include "specfem/setup.hpp" +#include "specfem/utilities/band.hpp" #include #include @@ -108,26 +108,20 @@ template struct AttenuationObjective { * type_real Q = 200.0; * type_real min_period = 0.01; * type_real max_period = 10.0; - * auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - * auto tau_eps = compute_tau_eps(Q, tau_sigma, min_frequency, - * max_frequency); + * specfem::utilities::Band period_band(min_period, + * max_period); auto tau_sigma = compute_tau_sigma(period_band); auto + * tau_eps = compute_tau_eps(Q, tau_sigma, period_band); * @endcode */ -template +template Kokkos::View compute_tau_eps( type_real Q, Kokkos::View tau_sigma, - type_real min_frequency, type_real max_frequency); + const specfem::utilities::Band band); } // namespace attenuation } // namespace specfem -extern template Kokkos::View - specfem::attenuation::compute_tau_eps( - type_real, - Kokkos::View, - type_real, type_real); +#include "compute_tau_eps.tpp" diff --git a/core/specfem/attenuation/compute_tau_eps.tpp b/core/specfem/attenuation/compute_tau_eps.tpp index 03b74e85f..d309bed71 100644 --- a/core/specfem/attenuation/compute_tau_eps.tpp +++ b/core/specfem/attenuation/compute_tau_eps.tpp @@ -1,10 +1,12 @@ #pragma once -#include "compute_tau_eps.hpp" #include "maxwell.hpp" +#include "compute_tau_eps.hpp" #include "specfem/constants.hpp" #include "specfem/optimization.hpp" #include "specfem/utilities/logspace.hpp" +#include "specfem/utilities/band.hpp" #include "specfem/setup.hpp" +#include "specfem/units.hpp" #include #include @@ -12,17 +14,23 @@ namespace specfem { namespace attenuation { -template +template Kokkos::View compute_tau_eps( type_real Q, Kokkos::View tau_sigma, - type_real min_frequency, type_real max_frequency) { + const specfem::utilities::Band band) { static_assert(N_SLS > 1, "N_SLS must be greater than 1 for tau_eps computation"); + // Convert to frequency band + const specfem::utilities::Band frequency_band(band); + + type_real min_frequency = frequency_band.min.raw(); + type_real max_frequency = frequency_band.max.raw(); + // Set up evaluation frequencies equally spaced in log10 const auto f = specfem::utilities::logspace(min_frequency, max_frequency); diff --git a/core/specfem/attenuation/compute_tau_sigma.hpp b/core/specfem/attenuation/compute_tau_sigma.hpp index e82e9a348..83ae8328e 100644 --- a/core/specfem/attenuation/compute_tau_sigma.hpp +++ b/core/specfem/attenuation/compute_tau_sigma.hpp @@ -2,6 +2,8 @@ #include "specfem/constants.hpp" #include "specfem/setup.hpp" +#include "specfem/units.hpp" +#include "specfem/utilities/band.hpp" #include #include @@ -32,14 +34,11 @@ namespace attenuation { * for (int i = 0; i < 4; ++i) std::cout << tau(i) << "\n"; * @endcode */ -template +template Kokkos::View -compute_tau_sigma(const type_real min_frequency, const type_real max_frequency); +compute_tau_sigma(const specfem::utilities::Band band); } // namespace attenuation } // namespace specfem -extern template Kokkos::View -specfem::attenuation::compute_tau_sigma( - const type_real, const type_real); +#include "compute_tau_sigma.tpp" diff --git a/core/specfem/attenuation/compute_tau_sigma.tpp b/core/specfem/attenuation/compute_tau_sigma.tpp index b3f42d3af..0d9fdc2af 100644 --- a/core/specfem/attenuation/compute_tau_sigma.tpp +++ b/core/specfem/attenuation/compute_tau_sigma.tpp @@ -1,20 +1,26 @@ #pragma once -#include "compute_tau_sigma.hpp" #include "specfem/constants.hpp" +#include "specfem/units.hpp" #include "specfem/setup.hpp" +#include "specfem/utilities/band.hpp" #include #include namespace specfem { namespace attenuation { -template +template Kokkos::View -compute_tau_sigma(const type_real min_frequency, const type_real max_frequency) { +compute_tau_sigma(const specfem::utilities::Band band) { static_assert(N_SLS > 1, "N_SLS must be greater than 1 to avoid division by zero"); - Kokkos::View tau_s( - "tau_sigma"); + // If the band is not convertible to a frequency band, this will throw an + // exception. This is intentional. + const specfem::utilities::Band frequency_band(band); + + // Extract the minimum and maximum frequencies from the band + type_real min_frequency = frequency_band.min.raw(); + type_real max_frequency = frequency_band.max.raw(); // logarithms of the input frequencies const type_real exp1 = std::log10(min_frequency); @@ -24,6 +30,9 @@ compute_tau_sigma(const type_real min_frequency, const type_real max_frequency) const type_real dexpval = (exp2 - exp1) / (static_cast(N_SLS) - 1); + Kokkos::View tau_s( + "tau_sigma"); + for (int i = 0; i < N_SLS; ++i) { tau_s(i) = 1.0 / (specfem::constants::pi * 2.0 * std::pow(10.0, exp1 + i * dexpval)); } diff --git a/core/specfem/data_access/check_compatibility.hpp b/core/specfem/data_access/check_compatibility.hpp index 743d301a5..b220af2ec 100644 --- a/core/specfem/data_access/check_compatibility.hpp +++ b/core/specfem/data_access/check_compatibility.hpp @@ -59,6 +59,20 @@ struct is_field_derivatives< specfem::data_access::DataClassType::field_derivatives> > : std::true_type {}; +/// @brief Detects the null/sentinel field-derivatives type +/// (specfem::point::null_field_derivatives). +template +struct is_null_field_derivatives : std::false_type {}; + +template +struct is_null_field_derivatives< + T, std::enable_if_t > : std::true_type {}; + +template +struct is_field_derivatives > + : std::true_type {}; + /// @brief Detects source term types template struct is_source : std::false_type {}; diff --git a/core/specfem/io/kernel/writer.hpp b/core/specfem/io/kernel/writer.hpp index d5a5105b2..86f6f50ed 100644 --- a/core/specfem/io/kernel/writer.hpp +++ b/core/specfem/io/kernel/writer.hpp @@ -1,7 +1,5 @@ #pragma once -#include "specfem/assembly/assembly.hpp" -#include "specfem/enums.hpp" #include "specfem/io/writer.hpp" namespace specfem { diff --git a/core/specfem/io/seismogram/reader.hpp b/core/specfem/io/seismogram/reader.hpp index f5bc32483..0fadba874 100644 --- a/core/specfem/io/seismogram/reader.hpp +++ b/core/specfem/io/seismogram/reader.hpp @@ -1,7 +1,8 @@ #pragma once #include "specfem/enums.hpp" -#include "specfem/io/reader.hpp" +#include "specfem/setup.hpp" +#include namespace specfem { namespace source_time_functions { diff --git a/core/specfem/io/wavefield/reader.tpp b/core/specfem/io/wavefield/reader.tpp index edc3f5d65..d4bd13c54 100644 --- a/core/specfem/io/wavefield/reader.tpp +++ b/core/specfem/io/wavefield/reader.tpp @@ -1,5 +1,6 @@ #pragma once +#include "specfem/assembly.hpp" #include "specfem/io/wavefield/reader.hpp" #include "specfem/utilities.hpp" #include "specfem/mpi.hpp" @@ -107,7 +108,6 @@ void specfem::io::wavefield_reader::run( specfem::assembly::assembly &assembly, const int istep) { auto &buffer = assembly.fields.buffer; - auto &boundary_values = assembly.boundary_values; typename IOLibrary::Group base_group = file.openGroup( std::string("/Step") + specfem::utilities::to_zero_lead(istep, 6)); diff --git a/core/specfem/medium/dim2/acoustic/isotropic/material.hpp b/core/specfem/medium/dim2/acoustic/isotropic/material.hpp index 61d1ce2b7..d25eb9266 100644 --- a/core/specfem/medium/dim2/acoustic/isotropic/material.hpp +++ b/core/specfem/medium/dim2/acoustic/isotropic/material.hpp @@ -146,9 +146,20 @@ class material > + template < + specfem::element::attenuation_tag T = AttenuationTag, + std::enable_if_t = 0> + inline specfem::point::properties< + specfem::tags::Tags > + get_properties() const { + return { static_cast(1.0) / static_cast(density), + this->kappa }; + } + + template < + specfem::element::attenuation_tag T = AttenuationTag, + std::enable_if_t< + T == specfem::element::attenuation_tag::constant_isotropic, int> = 0> inline specfem::point::properties< specfem::tags::Tags > get_properties() const { diff --git a/core/specfem/medium/dim2/elastic/isotropic/material.hpp b/core/specfem/medium/dim2/elastic/isotropic/material.hpp index 89349973b..8423cb4c4 100644 --- a/core/specfem/medium/dim2/elastic/isotropic/material.hpp +++ b/core/specfem/medium/dim2/elastic/isotropic/material.hpp @@ -1,8 +1,13 @@ #pragma once +#include "specfem/attenuation.hpp" +#include "specfem/constants.hpp" #include "specfem/element.hpp" #include "specfem/setup.hpp" #include "specfem/tags.hpp" +#include "specfem/units.hpp" +#include "specfem/utilities/band.hpp" +#include #include #include #include @@ -13,6 +18,61 @@ namespace medium_container { namespace impl { +template +struct ComputedAttenuationValues< + DimensionTag, MediumTag, + specfem::element::attenuation_tag::constant_isotropic, + std::enable_if_t::value> > { + + using view_type = Kokkos::View; + + type_real kappa_scale; ///< Scaling factor for bulk modulus attenuation + type_real mu_scale; ///< Scaling factor for shear modulus attenuation + + view_type tau_epsilon_kappa; ///< Relaxation times for SLS mechanisms + view_type tau_epsilon_mu; ///< Relaxation times for SLS mechanisms + specfem::attenuation::AttenuationPropertyValues + kappa_attenuation_properties; ///< Struct to hold computed attenuation + ///< properties + specfem::attenuation::AttenuationPropertyValues + mu_attenuation_properties; ///< Struct to hold computed attenuation + ///< properties + +public: + ComputedAttenuationValues() = default; + + ComputedAttenuationValues( + const type_real &Qkappa, const type_real &Qmu, const type_real &f0, + const type_real &fc, + const specfem::utilities::Band &band, + const view_type &tau_sigma) { + + this->tau_epsilon_kappa = + specfem::attenuation::compute_tau_eps( + Qkappa, tau_sigma, band); + + this->tau_epsilon_mu = + specfem::attenuation::compute_tau_eps( + Qmu, tau_sigma, band); + + this->kappa_attenuation_properties = + specfem::attenuation::get_attenuation_property_values< + specfem::constants::N_SLS>(tau_sigma, tau_epsilon_kappa); + + this->mu_attenuation_properties = + specfem::attenuation::get_attenuation_property_values< + specfem::constants::N_SLS>(tau_sigma, tau_epsilon_mu); + + this->kappa_scale = specfem::attenuation::get_attenuation_scale_factor< + specfem::constants::N_SLS>(fc, tau_epsilon_kappa, tau_sigma, Qkappa, + f0); + this->mu_scale = specfem::attenuation::get_attenuation_scale_factor< + specfem::constants::N_SLS>(fc, tau_epsilon_mu, tau_sigma, Qmu, f0); + } +}; + template struct AttenuationValues< @@ -23,6 +83,23 @@ struct AttenuationValues< type_real Qkappa; ///< Attenuation factor for bulk modulus type_real Qmu; ///< Attenuation factor for shear modulus + constexpr static auto dimension_tag = + DimensionTag; ///< Dimension of the material + constexpr static auto medium_tag = MediumTag; ///< Medium tag + constexpr static auto attenuation_tag = + specfem::element::attenuation_tag::constant_isotropic; ///< Attenuation + ///< tag +private: + ComputedAttenuationValues + computed_values; ///< Struct to hold computed attenuation properties + + bool compute_properties_called = + false; ///< Flag to check if properties have been computed + + using view_type = Kokkos::View; + +public: AttenuationValues() = default; AttenuationValues(const type_real &Qkappa, const type_real &Qmu) @@ -39,11 +116,50 @@ struct AttenuationValues< std::abs(this->Qmu - other.Qmu) < 1e-6); } + const ComputedAttenuationValues & + compute_attenuation_properties( + const type_real &f0, const type_real &fc, + const specfem::utilities::Band &band, + const view_type &tau_sigma) { + + if (this->compute_properties_called) { + return this->computed_values; + } + + this->computed_values = + ComputedAttenuationValues( + this->Qkappa, this->Qmu, f0, fc, band, tau_sigma); + + this->compute_properties_called = true; + + return this->computed_values; + } + + // Getters for computed properties with error handling + const ComputedAttenuationValues & + get_computed_values() const { + if (!this->compute_properties_called) { + throw std::runtime_error( + "Attenuation properties have not been computed yet. Call " + "compute_attenuation_properties() first."); + } + return this->computed_values; + } + std::string print() const { std::ostringstream message; message << " Qkappa : " << this->Qkappa << "\n" << " Qmu : " << this->Qmu << "\n"; + if (this->compute_properties_called) { + message << " Computed attenuation properties:\n" + << " kappa scale factor: " + << this->computed_values.kappa_scale << "\n" + << " mu scale factor: " + << this->computed_values.mu_scale << "\n"; + } else { + message << " Attenuation properties have not been computed yet.\n"; + } return message.str(); } @@ -178,12 +294,33 @@ struct material< * * @return specfem::point::properties Material properties */ + template < + specfem::element::attenuation_tag T = AttenuationTag, + std::enable_if_t = 0> inline specfem::point::properties< specfem::tags::Tags > get_properties() const { return { this->kappa, this->mu, this->density }; } + /** + * @brief Get the material properties with attenuation scaling + * + * @return specfem::point::properties Material properties with attenuation + * scaling + */ + template < + specfem::element::attenuation_tag T = AttenuationTag, + std::enable_if_t< + T == specfem::element::attenuation_tag::constant_isotropic, int> = 0> + inline specfem::point::properties< + specfem::tags::Tags > + get_properties() const { + return { this->kappa * attenuation::get_computed_values().kappa_scale, + this->mu * attenuation::get_computed_values().mu_scale, + this->density }; + } + /** * @brief Print the material properties * diff --git a/core/specfem/medium_container/impl/attenuation_values.hpp b/core/specfem/medium_container/impl/attenuation_values.hpp index 9e4d4f6c2..07bb66950 100644 --- a/core/specfem/medium_container/impl/attenuation_values.hpp +++ b/core/specfem/medium_container/impl/attenuation_values.hpp @@ -21,4 +21,11 @@ class AttenuationValues +class ComputedAttenuationValues; + } // namespace specfem::medium_container::impl diff --git a/core/specfem/medium_container/material.hpp b/core/specfem/medium_container/material.hpp index c62489988..7e5effccf 100644 --- a/core/specfem/medium_container/material.hpp +++ b/core/specfem/medium_container/material.hpp @@ -58,6 +58,5 @@ class material; #include "specfem/medium/dim2/elastic/isotropic_cosserat/material.hpp" #include "specfem/medium/dim2/electromagnetic/isotropic/material.hpp" #include "specfem/medium/dim2/poroelastic/isotropic/material.hpp" - #include "specfem/medium/dim3/acoustic/isotropic/material.hpp" #include "specfem/medium/dim3/elastic/isotropic/material.hpp" diff --git a/core/specfem/point/attenuation.hpp b/core/specfem/point/attenuation.hpp index 60a34d6bf..c69dbf6cd 100644 --- a/core/specfem/point/attenuation.hpp +++ b/core/specfem/point/attenuation.hpp @@ -108,6 +108,10 @@ struct attenuationkappa_common_factor; - result.mu_common_factor = this->mu_common_factor; + result.kappa_relaxation_rate = this->kappa_relaxation_rate; + result.mu_relaxation_rate = this->mu_relaxation_rate; result.alpha_rk = this->alpha_rk; result.beta_rk = this->beta_rk; result.gamma_rk = this->gamma_rk; @@ -213,8 +217,8 @@ struct attenuationkappa_common_factor; - result.mu_common_factor = this->mu_common_factor; + result.kappa_relaxation_rate = this->kappa_relaxation_rate; + result.mu_relaxation_rate = this->mu_relaxation_rate; result.alpha_rk = this->alpha_rk; result.beta_rk = this->beta_rk; result.gamma_rk = this->gamma_rk; @@ -231,9 +235,9 @@ struct attenuationkappa_common_factor; - result.mu_common_factor = this->mu_common_factor; + result.kappa_relaxation_rate = this->kappa_relaxation_rate; + result.mu_relaxation_rate = this->mu_relaxation_rate; result.alpha_rk = this->alpha_rk; result.beta_rk = this->beta_rk; result.gamma_rk = this->gamma_rk; @@ -434,8 +442,8 @@ struct attenuationkappa_common_factor; - result.mu_common_factor = this->mu_common_factor; + result.kappa_relaxation_rate = this->kappa_relaxation_rate; + result.mu_relaxation_rate = this->mu_relaxation_rate; result.alpha_rk = this->alpha_rk; result.beta_rk = this->beta_rk; result.gamma_rk = this->gamma_rk; @@ -455,9 +463,9 @@ struct attenuation #include @@ -48,15 +50,18 @@ bool specfem::program::qplots(const type_real Q, const type_real minfreq, specfem::Logger::info("Directory already exists: " + qplot_dir.string()); } + using namespace specfem::units::unit_symbols; + const specfem::utilities::Band band(minfreq * Hz, + maxfreq * Hz); + // Get tau_sigma const auto tau_sigma = - specfem::attenuation::compute_tau_sigma( - minfreq, maxfreq); + specfem::attenuation::compute_tau_sigma(band); // Compute tau_eps for the given Q const auto tau_eps = specfem::attenuation::compute_tau_eps( - Q, tau_sigma, minfreq, maxfreq); + Q, tau_sigma, band); // Print tau_sigma and tau_eps for debugging specfem::Logger::info("Computed tau_sigma:"); diff --git a/core/specfem/source_time_functions/external.cpp b/core/specfem/source_time_functions/external.cpp index 22bf95e7f..0675e3dcc 100644 --- a/core/specfem/source_time_functions/external.cpp +++ b/core/specfem/source_time_functions/external.cpp @@ -1,7 +1,7 @@ #include "external.hpp" #include "specfem/enums.hpp" -#include "specfem/io.hpp" +#include "specfem/io/seismogram/reader.hpp" #include "specfem/utilities.hpp" #include #include diff --git a/core/specfem/units.hpp b/core/specfem/units.hpp index 57e7e64a9..c1ed7a372 100644 --- a/core/specfem/units.hpp +++ b/core/specfem/units.hpp @@ -3,8 +3,6 @@ #include "units/parse.hpp" #include "units/quantity.hpp" -#include - /** * @namespace specfem::units * @brief Type-safe physical units and conversions. @@ -32,15 +30,15 @@ namespace specfem::units::unit_symbols { /// Tag for constructing Meters struct meter_tag {}; -KOKKOS_FORCEINLINE_FUNCTION constexpr meter_tag m{}; +constexpr meter_tag m{}; /// Tag for constructing Kilometers struct kilometer_tag {}; -KOKKOS_FORCEINLINE_FUNCTION constexpr kilometer_tag km{}; +constexpr kilometer_tag km{}; /// Tag for constructing Seconds struct second_tag {}; -KOKKOS_FORCEINLINE_FUNCTION constexpr second_tag s{}; +constexpr second_tag s{}; template constexpr Meters operator*(N v, meter_tag) { return Meters(type_real(v)); @@ -63,11 +61,11 @@ template constexpr Seconds operator*(second_tag, N v) { /// Tag for constructing MetersPerSecond struct m_per_s_tag {}; -KOKKOS_FORCEINLINE_FUNCTION constexpr m_per_s_tag mps{}; +constexpr m_per_s_tag mps{}; /// Tag for constructing KilometersPerSecond struct km_per_s_tag {}; -KOKKOS_FORCEINLINE_FUNCTION constexpr km_per_s_tag kmps{}; +constexpr km_per_s_tag kmps{}; template constexpr MetersPerSecond operator*(N v, m_per_s_tag) { return MetersPerSecond(type_real(v)); @@ -84,4 +82,15 @@ constexpr KilometersPerSecond operator*(km_per_s_tag, N v) { return KilometersPerSecond(type_real(v)); } +// Tag for constructing Hertz +struct hertz_tag {}; +constexpr hertz_tag Hz{}; + +template constexpr Hertz operator*(N v, hertz_tag) { + return Hertz(type_real(v)); +} +template constexpr Hertz operator*(hertz_tag, N v) { + return Hertz(type_real(v)); +} + } // namespace specfem::units::unit_symbols diff --git a/core/specfem/units/conversions.hpp b/core/specfem/units/conversions.hpp index 4c62b3888..9b0b7c03e 100644 --- a/core/specfem/units/conversions.hpp +++ b/core/specfem/units/conversions.hpp @@ -1,5 +1,6 @@ #pragma once #include "quantity.hpp" +#include #include namespace specfem::units { @@ -37,14 +38,15 @@ template struct unit_cast_impl { /// Identity conversion (no-op) template struct unit_cast_impl { - static constexpr T call(T q) noexcept { return q; } + KOKKOS_FUNCTION static constexpr T call(T q) noexcept { return q; } }; /// Scale conversion within same dimension (e.g., meters ↔ kilometers) template struct unit_cast_impl, Quantity, std::enable_if_t::value> > { - static constexpr Quantity call(Quantity q) noexcept { + KOKKOS_FUNCTION static constexpr Quantity + call(Quantity q) noexcept { constexpr type_real factor = ratio_value >; return Quantity(q.raw() * factor); } @@ -63,42 +65,42 @@ struct unit_cast_impl, Quantity, /// Seconds → Hertz template <> struct unit_cast_impl { - static constexpr Hertz call(Seconds s) noexcept { + KOKKOS_FUNCTION static constexpr Hertz call(Seconds s) noexcept { return Hertz(type_real(1) / s.raw()); } }; /// Hertz → Seconds template <> struct unit_cast_impl { - static constexpr Seconds call(Hertz f) noexcept { + KOKKOS_FUNCTION static constexpr Seconds call(Hertz f) noexcept { return Seconds(type_real(1) / f.raw()); } }; /// Seconds → Omega template <> struct unit_cast_impl { - static constexpr Omega call(Seconds s) noexcept { + KOKKOS_FUNCTION static constexpr Omega call(Seconds s) noexcept { return Omega(impl::two_pi / s.raw()); } }; /// Omega → Seconds template <> struct unit_cast_impl { - static constexpr Seconds call(Omega w) noexcept { + KOKKOS_FUNCTION static constexpr Seconds call(Omega w) noexcept { return Seconds(impl::two_pi / w.raw()); } }; /// Hertz → Omega template <> struct unit_cast_impl { - static constexpr Omega call(Hertz f) noexcept { + KOKKOS_FUNCTION static constexpr Omega call(Hertz f) noexcept { return Omega(f.raw() * impl::two_pi); } }; /// Omega → Hertz template <> struct unit_cast_impl { - static constexpr Hertz call(Omega w) noexcept { + KOKKOS_FUNCTION static constexpr Hertz call(Omega w) noexcept { return Hertz(w.raw() / impl::two_pi); } }; @@ -125,7 +127,8 @@ template <> struct unit_cast_impl { * auto freq = unit_cast(Seconds(0.5)); // 2.0 Hz * @endcode */ -template constexpr To unit_cast(From q) { +template +KOKKOS_FUNCTION constexpr To unit_cast(From q) { return unit_cast_impl::call(q); } diff --git a/core/specfem/units/quantity.hpp b/core/specfem/units/quantity.hpp index 761546dbc..9e2508568 100644 --- a/core/specfem/units/quantity.hpp +++ b/core/specfem/units/quantity.hpp @@ -1,10 +1,9 @@ #pragma once #include "specfem/setup.hpp" +#include #include #include -#include - namespace specfem::units { /** @@ -68,7 +67,7 @@ using ratio_gcd = * @return constexpr type_real Ratio as a floating-point value */ template -KOKKOS_FORCEINLINE_FUNCTION constexpr type_real ratio_value = +constexpr type_real ratio_value = static_cast(R::num) / static_cast(R::den); /** @@ -96,7 +95,8 @@ template > class Quantity { * * @param v Value in the specified units (default: 0.0) */ - constexpr explicit Quantity(type_real v = 0.0) noexcept : value_(v) {} + KOKKOS_FUNCTION constexpr explicit Quantity(type_real v = 0.0) noexcept + : value_(v) {} /** * @brief Extract the raw numeric value. @@ -105,60 +105,75 @@ template > class Quantity { * * @return constexpr type_real Numeric value without units */ - [[nodiscard]] constexpr type_real raw() const noexcept { return value_; } + [[nodiscard]] KOKKOS_FUNCTION constexpr type_real raw() const noexcept { + return value_; + } /** * @brief Explicit conversion to raw numeric value. * * @return constexpr type_real Numeric value without units */ - explicit constexpr operator type_real() const noexcept { return value_; } + KOKKOS_FUNCTION explicit constexpr operator type_real() const noexcept { + return value_; + } /// Unary negation - constexpr Quantity operator-() const noexcept { return Quantity(-value_); } + KOKKOS_FUNCTION constexpr Quantity operator-() const noexcept { + return Quantity(-value_); + } // Arithmetic within the same dimension and scale - constexpr Quantity operator+(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr Quantity operator+(Quantity o) const noexcept { return Quantity(value_ + o.value_); } - constexpr Quantity operator-(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr Quantity operator-(Quantity o) const noexcept { return Quantity(value_ - o.value_); } - constexpr Quantity operator*(type_real s) const noexcept { + KOKKOS_FUNCTION constexpr Quantity operator*(type_real s) const noexcept { return Quantity(value_ * s); } - constexpr Quantity operator/(type_real s) const noexcept { + KOKKOS_FUNCTION constexpr Quantity operator/(type_real s) const noexcept { return Quantity(value_ / s); } // Comparisons (same scale) - constexpr bool operator==(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr bool operator==(Quantity o) const noexcept { return value_ == o.value_; } - constexpr bool operator!=(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr bool operator!=(Quantity o) const noexcept { return value_ != o.value_; } - constexpr bool operator<(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr bool operator<(Quantity o) const noexcept { return value_ < o.value_; } - constexpr bool operator<=(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr bool operator<=(Quantity o) const noexcept { return value_ <= o.value_; } - constexpr bool operator>(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr bool operator>(Quantity o) const noexcept { return value_ > o.value_; } - constexpr bool operator>=(Quantity o) const noexcept { + KOKKOS_FUNCTION constexpr bool operator>=(Quantity o) const noexcept { return value_ >= o.value_; } }; /// Scalar multiplication from left template -constexpr Quantity operator*(type_real s, - Quantity q) noexcept { +KOKKOS_FUNCTION constexpr Quantity +operator*(type_real s, Quantity q) noexcept { return q * s; } +/// Scalar division from left: s / Quantity,S> -> +/// Quantity,S> +template +KOKKOS_FUNCTION constexpr auto operator/(type_real s, + Quantity, S> q) + -> Quantity, S> { + return Quantity, S>(s / q.raw()); +} + /** * @name Mixed-scale arithmetic * @brief Operations between quantities with different scales but same @@ -170,7 +185,7 @@ constexpr Quantity operator*(type_real s, template ::value> > -constexpr auto operator+(Quantity a, Quantity b) +KOKKOS_FUNCTION constexpr auto operator+(Quantity a, Quantity b) -> Quantity > { using Rgcd = impl::ratio_gcd; return Quantity(a.raw() * ratio_value > + @@ -179,7 +194,7 @@ constexpr auto operator+(Quantity a, Quantity b) template ::value> > -constexpr auto operator-(Quantity a, Quantity b) +KOKKOS_FUNCTION constexpr auto operator-(Quantity a, Quantity b) -> Quantity > { using Rgcd = impl::ratio_gcd; return Quantity(a.raw() * ratio_value > - @@ -197,7 +212,8 @@ constexpr auto operator-(Quantity a, Quantity b) template ::value> > -constexpr bool operator==(Quantity a, Quantity b) noexcept { +KOKKOS_FUNCTION constexpr bool operator==(Quantity a, + Quantity b) noexcept { using Rgcd = impl::ratio_gcd; return a.raw() * ratio_value > == b.raw() * ratio_value >; @@ -205,13 +221,15 @@ constexpr bool operator==(Quantity a, Quantity b) noexcept { template ::value> > -constexpr bool operator!=(Quantity a, Quantity b) noexcept { +KOKKOS_FUNCTION constexpr bool operator!=(Quantity a, + Quantity b) noexcept { return !(a == b); } template ::value> > -constexpr bool operator<(Quantity a, Quantity b) noexcept { +KOKKOS_FUNCTION constexpr bool operator<(Quantity a, + Quantity b) noexcept { using Rgcd = impl::ratio_gcd; return a.raw() * ratio_value > < b.raw() * ratio_value >; @@ -219,19 +237,22 @@ constexpr bool operator<(Quantity a, Quantity b) noexcept { template ::value> > -constexpr bool operator<=(Quantity a, Quantity b) noexcept { +KOKKOS_FUNCTION constexpr bool operator<=(Quantity a, + Quantity b) noexcept { return !(b < a); } template ::value> > -constexpr bool operator>(Quantity a, Quantity b) noexcept { +KOKKOS_FUNCTION constexpr bool operator>(Quantity a, + Quantity b) noexcept { return b < a; } template ::value> > -constexpr bool operator>=(Quantity a, Quantity b) noexcept { +KOKKOS_FUNCTION constexpr bool operator>=(Quantity a, + Quantity b) noexcept { return !(a < b); } @@ -248,8 +269,8 @@ constexpr bool operator>=(Quantity a, Quantity b) noexcept { template -constexpr auto operator*(Quantity, S1> a, - Quantity, S2> b) +KOKKOS_FUNCTION constexpr auto operator*(Quantity, S1> a, + Quantity, S2> b) -> Quantity, std::ratio_multiply > { return Quantity, @@ -258,8 +279,8 @@ constexpr auto operator*(Quantity, S1> a, template -constexpr auto operator/(Quantity, S1> a, - Quantity, S2> b) +KOKKOS_FUNCTION constexpr auto operator/(Quantity, S1> a, + Quantity, S2> b) -> Quantity, std::ratio_divide > { return Quantity, diff --git a/tests/unit-tests/attenuation/compute_band_tests.cpp b/tests/unit-tests/attenuation/compute_band_tests.cpp index 9b5b2b9f6..b69162a2f 100644 --- a/tests/unit-tests/attenuation/compute_band_tests.cpp +++ b/tests/unit-tests/attenuation/compute_band_tests.cpp @@ -6,6 +6,7 @@ #include using specfem::attenuation::compute_band; +using specfem::units::Hertz; using specfem::units::Omega; using specfem::units::Seconds; using specfem::units::unit_cast; @@ -16,52 +17,42 @@ using specfem::utilities::is_close; // compute_band tests // --------------------------------------------------------------------------- -// min_period is always equal to the input min_resolved_period -TEST(Attenuation_ComputeBand, MinPeriodEqualsInput) { +// max frequency equals 1/min_resolved_period +TEST(Attenuation_ComputeBand, MaxFreqEqualsInverseInput) { constexpr type_real min_resolved = 0.5; - EXPECT_TRUE(is_close( - unit_cast(compute_band<2>(Seconds(min_resolved))).min.raw(), - min_resolved)); - EXPECT_TRUE(is_close( - unit_cast(compute_band<3>(Seconds(min_resolved))).min.raw(), - min_resolved)); - EXPECT_TRUE(is_close( - unit_cast(compute_band<4>(Seconds(min_resolved))).min.raw(), - min_resolved)); - EXPECT_TRUE(is_close( - unit_cast(compute_band<5>(Seconds(min_resolved))).min.raw(), - min_resolved)); + EXPECT_TRUE(is_close(compute_band<2>(Seconds(min_resolved)).max.raw(), + 1.0 / min_resolved)); + EXPECT_TRUE(is_close(compute_band<3>(Seconds(min_resolved)).max.raw(), + 1.0 / min_resolved)); + EXPECT_TRUE(is_close(compute_band<4>(Seconds(min_resolved)).max.raw(), + 1.0 / min_resolved)); + EXPECT_TRUE(is_close(compute_band<5>(Seconds(min_resolved)).max.raw(), + 1.0 / min_resolved)); } -// max_period must be strictly larger than min_period for all valid N_SLS -TEST(Attenuation_ComputeBand, MaxPeriodExceedsMin) { +// max frequency must be strictly larger than min frequency for all valid N_SLS +TEST(Attenuation_ComputeBand, MaxFreqExceedsMin) { constexpr type_real min_resolved = 0.5; - EXPECT_GT( - unit_cast(compute_band<2>(Seconds(min_resolved))).max.raw(), - min_resolved); - EXPECT_GT( - unit_cast(compute_band<3>(Seconds(min_resolved))).max.raw(), - min_resolved); - EXPECT_GT( - unit_cast(compute_band<4>(Seconds(min_resolved))).max.raw(), - min_resolved); - EXPECT_GT( - unit_cast(compute_band<5>(Seconds(min_resolved))).max.raw(), - min_resolved); + EXPECT_GT(compute_band<2>(Seconds(min_resolved)).max.raw(), + compute_band<2>(Seconds(min_resolved)).min.raw()); + EXPECT_GT(compute_band<3>(Seconds(min_resolved)).max.raw(), + compute_band<3>(Seconds(min_resolved)).min.raw()); + EXPECT_GT(compute_band<4>(Seconds(min_resolved)).max.raw(), + compute_band<4>(Seconds(min_resolved)).min.raw()); + EXPECT_GT(compute_band<5>(Seconds(min_resolved)).max.raw(), + compute_band<5>(Seconds(min_resolved)).min.raw()); } -// The decade width theta(N_SLS) gives max_period = min_period * 10^theta +// The decade width theta(N_SLS) gives max_freq/min_freq = 10^theta // Verified values: theta = {0.75, 1.75, 2.25, 2.85} for N_SLS = {2,3,4,5} TEST(Attenuation_ComputeBand, DecadeWidthMatchesTheta) { constexpr type_real min_resolved = 1.0; auto check = [&](auto band, double theta) { - auto period_band = unit_cast(band); type_real expected_ratio = static_cast(std::pow(10.0, theta)); - EXPECT_TRUE( - is_close((period_band.max / period_band.min).raw(), expected_ratio)); + EXPECT_TRUE(is_close((band.max / band.min).raw(), expected_ratio)); }; check(compute_band<2>(Seconds(min_resolved)), 0.75); @@ -70,28 +61,28 @@ TEST(Attenuation_ComputeBand, DecadeWidthMatchesTheta) { check(compute_band<5>(Seconds(min_resolved)), 2.85); } -// The band scales linearly with min_resolved_period +// The band scales inversely with min_resolved_period TEST(Attenuation_ComputeBand, ScalesWithInput) { constexpr type_real base = 0.1; constexpr type_real scale = 5.0; - auto b1 = unit_cast(compute_band<3>(Seconds(base))); - auto b2 = unit_cast(compute_band<3>(Seconds(base * scale))); + auto b1 = compute_band<3>(Seconds(base)); + auto b2 = compute_band<3>(Seconds(base * scale)); - EXPECT_TRUE(is_close(b2.min.raw(), b1.min.raw() * scale)); - EXPECT_TRUE(is_close(b2.max.raw(), b1.max.raw() * scale)); + EXPECT_TRUE(is_close(b2.min.raw(), b1.min.raw() / scale)); + EXPECT_TRUE(is_close(b2.max.raw(), b1.max.raw() / scale)); } -// Increasing N_SLS widens the band (larger theta -> larger max_period) +// Increasing N_SLS widens the band (larger theta -> lower min frequency) TEST(Attenuation_ComputeBand, MoreSLSGivesWiderBand) { constexpr type_real min_resolved = 1.0; - auto b2 = unit_cast(compute_band<2>(Seconds(min_resolved))); - auto b3 = unit_cast(compute_band<3>(Seconds(min_resolved))); - auto b4 = unit_cast(compute_band<4>(Seconds(min_resolved))); - auto b5 = unit_cast(compute_band<5>(Seconds(min_resolved))); + auto b2 = compute_band<2>(Seconds(min_resolved)); + auto b3 = compute_band<3>(Seconds(min_resolved)); + auto b4 = compute_band<4>(Seconds(min_resolved)); + auto b5 = compute_band<5>(Seconds(min_resolved)); - EXPECT_LT(b2.max, b3.max); - EXPECT_LT(b3.max, b4.max); - EXPECT_LT(b4.max, b5.max); + EXPECT_GT(b2.min, b3.min); + EXPECT_GT(b3.min, b4.min); + EXPECT_GT(b4.min, b5.min); } diff --git a/tests/unit-tests/attenuation/compute_tau_eps_tests.cpp b/tests/unit-tests/attenuation/compute_tau_eps_tests.cpp index 201852080..d79e30853 100644 --- a/tests/unit-tests/attenuation/compute_tau_eps_tests.cpp +++ b/tests/unit-tests/attenuation/compute_tau_eps_tests.cpp @@ -2,6 +2,8 @@ #include "specfem/attenuation/compute_tau_eps.tpp" #include "specfem/attenuation/compute_tau_sigma.tpp" #include "specfem/constants.hpp" +#include "specfem/units.hpp" +#include "specfem/utilities/band.hpp" #include #include @@ -9,6 +11,7 @@ using specfem::attenuation::compute_tau_eps; using specfem::attenuation::compute_tau_sigma; using specfem::attenuation::maxwell; using specfem::constants::NF_ATTENUATION; +using namespace specfem::units::unit_symbols; // Helper function to compute achieved Q from tau_eps and tau_sigma template @@ -137,10 +140,11 @@ TEST(Attenuation_ComputeTauEps, Savage2010_NSLS3_1p7Decades) { // 1.7 decades: max_frequency=10 Hz, min_frequency=10/10^1.7 const type_real max_frequency = 10.0; const type_real min_frequency = max_frequency / std::pow(10.0, 1.7); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error(tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); @@ -168,10 +172,11 @@ TEST(Attenuation_ComputeTauEps, Savage2010_NSLS2_0p9Decades) { // 0.9 decades: max_frequency=10 Hz, min_frequency=10/10^0.9 const type_real max_frequency = 10.0; const type_real min_frequency = max_frequency / std::pow(10.0, 0.9); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error(tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); @@ -193,10 +198,11 @@ TEST(Attenuation_ComputeTauEps, Savage2010_NSLS4_2p5Decades) { // 2.5 decades: max_frequency=10 Hz, min_frequency=10/10^2.5 const type_real max_frequency = 10.0; const type_real min_frequency = max_frequency / std::pow(10.0, 2.5); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error(tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); @@ -218,10 +224,11 @@ TEST(Attenuation_ComputeTauEps, Savage2010_NSLS5_3p0Decades) { // 3.0 decades: max_frequency=10 Hz, min_frequency=10/10^3.0=0.01 Hz const type_real max_frequency = 10.0; const type_real min_frequency = max_frequency / std::pow(10.0, 3.0); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error(tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); @@ -244,9 +251,12 @@ TEST(Attenuation_ComputeTauEps, Savage2010_ErrorIncreasesWithBandwidth) { // Compute error at 1.7 decades (optimal for N_SLS=3) type_real min_frequency_optimal = max_frequency / std::pow(10.0, 1.7); auto tau_sigma_opt = - compute_tau_sigma(min_frequency_optimal, max_frequency); + compute_tau_sigma(specfem::utilities::Band( + min_frequency_optimal * Hz, max_frequency * Hz)); auto tau_eps_opt = compute_tau_eps( - target_Q, tau_sigma_opt, min_frequency_optimal, max_frequency); + target_Q, tau_sigma_opt, + specfem::utilities::Band( + min_frequency_optimal * Hz, max_frequency * Hz)); type_real error_optimal = compute_lsq_error(tau_sigma_opt, tau_eps_opt, target_Q, min_frequency_optimal, max_frequency); @@ -254,9 +264,12 @@ TEST(Attenuation_ComputeTauEps, Savage2010_ErrorIncreasesWithBandwidth) { // Compute error at 3.0 decades (too wide for N_SLS=3) type_real min_frequency_wide = max_frequency / std::pow(10.0, 3.0); auto tau_sigma_wide = - compute_tau_sigma(min_frequency_wide, max_frequency); - auto tau_eps_wide = compute_tau_eps(target_Q, tau_sigma_wide, - min_frequency_wide, max_frequency); + compute_tau_sigma(specfem::utilities::Band( + min_frequency_wide * Hz, max_frequency * Hz)); + auto tau_eps_wide = + compute_tau_eps(target_Q, tau_sigma_wide, + specfem::utilities::Band( + min_frequency_wide * Hz, max_frequency * Hz)); type_real error_wide = compute_lsq_error(tau_sigma_wide, tau_eps_wide, target_Q, min_frequency_wide, max_frequency); @@ -276,21 +289,21 @@ TEST(Attenuation_ComputeTauEps, Savage2010_MoreSLSReducesError) { // 2 decades constexpr type_real max_frequency = 10.0; constexpr type_real min_frequency = 0.1; // 2 decades + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); // N_SLS = 2 { constexpr int N_SLS = 2; - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real error_2 = compute_lsq_error(tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); // N_SLS = 4 constexpr int N_SLS_4 = 4; - auto tau_sigma_4 = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps_4 = compute_tau_eps(target_Q, tau_sigma_4, - min_frequency, max_frequency); + auto tau_sigma_4 = compute_tau_sigma(band); + auto tau_eps_4 = compute_tau_eps(target_Q, tau_sigma_4, band); type_real error_4 = compute_lsq_error( tau_sigma_4, tau_eps_4, target_Q, min_frequency, max_frequency); @@ -312,10 +325,11 @@ TEST(Attenuation_ComputeTauEps, AchievesTargetQ_Medium) { constexpr type_real target_Q = 200.0; constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real achieved_Q = compute_achieved_Q( tau_sigma, tau_eps, min_frequency, max_frequency); @@ -331,10 +345,11 @@ TEST(Attenuation_ComputeTauEps, AchievesTargetQ_Low) { constexpr type_real target_Q = 50.0; constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real achieved_Q = compute_achieved_Q( tau_sigma, tau_eps, min_frequency, max_frequency); @@ -349,10 +364,11 @@ TEST(Attenuation_ComputeTauEps, AchievesTargetQ_High) { constexpr type_real target_Q = 1000.0; constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real achieved_Q = compute_achieved_Q( tau_sigma, tau_eps, min_frequency, max_frequency); @@ -367,10 +383,11 @@ TEST(Attenuation_ComputeTauEps, TauEpsGreaterThanTauS) { constexpr type_real target_Q = 200.0; constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); for (int j = 0; j < N_SLS; ++j) { EXPECT_GT(tau_eps(j), tau_sigma(j)) @@ -384,10 +401,11 @@ TEST(Attenuation_ComputeTauEps, ReturnsCorrectSize) { constexpr type_real target_Q = 200.0; constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); EXPECT_EQ(tau_eps.extent(0), N_SLS); } @@ -400,6 +418,8 @@ TEST(Attenuation_ComputeTauEps, DifferentNumberOfSLS) { // 3 decades: min_frequency=0.1, max_frequency=100.0 constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); type_real bandwidth = bandwidth_decades(min_frequency, max_frequency); std::cout << "Testing N_SLS comparison at " << bandwidth << " decades" @@ -410,9 +430,8 @@ TEST(Attenuation_ComputeTauEps, DifferentNumberOfSLS) { // At 3 decades, expect much higher error { constexpr int N_SLS = 2; - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error( tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); // At 3 decades, N_SLS=2 should have significant error (>10% per Savage @@ -428,9 +447,8 @@ TEST(Attenuation_ComputeTauEps, DifferentNumberOfSLS) { // At 3 decades, expect slightly higher but still reasonable error { constexpr int N_SLS = 4; - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error( tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); EXPECT_LT(lsq_error * 100.0, 5.0) @@ -443,9 +461,8 @@ TEST(Attenuation_ComputeTauEps, DifferentNumberOfSLS) { // Savage 2010: N_SLS=5 achieves ~1% error at ~3.0 decades { constexpr int N_SLS = 5; - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error( tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); EXPECT_LT(lsq_error * 100.0, 2.0) @@ -467,9 +484,10 @@ TEST(Attenuation_ComputeTauEps, DifferentFrequencyRanges) { { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 1.0; // 1 decade - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error( tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); // 1 decade is well below optimal 1.7 decades, should have very low error @@ -483,9 +501,10 @@ TEST(Attenuation_ComputeTauEps, DifferentFrequencyRanges) { { constexpr type_real min_frequency = 1.0 / 5.01; // ~1.7 decades constexpr type_real max_frequency = 10.0; - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error( tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); EXPECT_LT(lsq_error * 100.0, 2.0) @@ -500,9 +519,10 @@ TEST(Attenuation_ComputeTauEps, DifferentFrequencyRanges) { { constexpr type_real min_frequency = 0.01; constexpr type_real max_frequency = 1000.0; // 5 decades - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = compute_tau_eps(target_Q, tau_sigma, min_frequency, - max_frequency); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); type_real lsq_error = compute_lsq_error( tau_sigma, tau_eps, target_Q, min_frequency, max_frequency); // 5 decades far exceeds optimal, expect significant error @@ -521,10 +541,11 @@ TEST(Attenuation_ComputeTauEps, QReconstructionOverFrequencies) { constexpr type_real target_Q = 200.0; constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; // 3 decades + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); // Set up evaluation frequencies Kokkos::View band(min_frequency * Hz, + max_frequency * Hz); - auto tau_sigma = compute_tau_sigma(min_frequency, max_frequency); - auto tau_eps = - compute_tau_eps(target_Q, tau_sigma, min_frequency, max_frequency); + auto tau_sigma = compute_tau_sigma(band); + auto tau_eps = compute_tau_eps(target_Q, tau_sigma, band); for (int j = 0; j < N_SLS; ++j) { EXPECT_TRUE(std::isfinite(tau_eps(j))) diff --git a/tests/unit-tests/attenuation/compute_tau_sigma_tests.cpp b/tests/unit-tests/attenuation/compute_tau_sigma_tests.cpp index 5343e0f7b..58cdb0d9e 100644 --- a/tests/unit-tests/attenuation/compute_tau_sigma_tests.cpp +++ b/tests/unit-tests/attenuation/compute_tau_sigma_tests.cpp @@ -1,19 +1,24 @@ #include "specfem/attenuation.hpp" #include "specfem/attenuation/compute_tau_sigma.tpp" +#include "specfem/units.hpp" +#include "specfem/utilities/band.hpp" #include #include using specfem::attenuation::compute_tau_sigma; +using namespace specfem::units::unit_symbols; // Test that the function returns the correct number of elements TEST(Attenuation_ComputeTauSigma, ReturnsCorrectSize) { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; - auto tau_s_3 = compute_tau_sigma<3>(min_frequency, max_frequency); + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); + auto tau_s_3 = compute_tau_sigma<3>(band); EXPECT_EQ(tau_s_3.extent(0), 3); - auto tau_s_5 = compute_tau_sigma<5>(min_frequency, max_frequency); + auto tau_s_5 = compute_tau_sigma<5>(band); EXPECT_EQ(tau_s_5.extent(0), 5); } @@ -21,8 +26,10 @@ TEST(Attenuation_ComputeTauSigma, ReturnsCorrectSize) { TEST(Attenuation_ComputeTauSigma, ValuesArePositive) { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_s = compute_tau_sigma<3>(min_frequency, max_frequency); + auto tau_s = compute_tau_sigma<3>(band); for (int i = 0; i < 3; ++i) { EXPECT_GT(tau_s(i), 0.0) << "tau_s(" << i << ") should be positive"; @@ -34,8 +41,10 @@ TEST(Attenuation_ComputeTauSigma, ValuesArePositive) { TEST(Attenuation_ComputeTauSigma, ValuesAreMonotonicallyDecreasing) { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_s = compute_tau_sigma<5>(min_frequency, max_frequency); + auto tau_s = compute_tau_sigma<5>(band); for (int i = 1; i < 5; ++i) { EXPECT_LT(tau_s(i), tau_s(i - 1)) @@ -48,8 +57,10 @@ TEST(Attenuation_ComputeTauSigma, EquallySpacedInLog10Frequency) { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; constexpr int N_SLS = 5; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_s = compute_tau_sigma(min_frequency, max_frequency); + auto tau_s = compute_tau_sigma(band); // Convert tau_s to frequencies: f = 1 / (2 * pi * tau_s) // Then check that log10(f) values are equally spaced @@ -73,8 +84,10 @@ TEST(Attenuation_ComputeTauSigma, BoundaryFrequenciesMatch) { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; constexpr int N_SLS = 3; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_s = compute_tau_sigma(min_frequency, max_frequency); + auto tau_s = compute_tau_sigma(band); // Convert tau_s to frequency: f = 1 / (2 * pi * tau_s) type_real freq_first = 1.0 / (2.0 * specfem::constants::pi * tau_s(0)); @@ -91,7 +104,8 @@ TEST(Attenuation_ComputeTauSigma, BoundaryFrequenciesMatch) { TEST(Attenuation_ComputeTauSigma, DifferentFrequencyRanges) { // Narrow range { - auto tau_s = compute_tau_sigma<3>(0.1, 1.0); + auto tau_s = compute_tau_sigma<3>( + specfem::utilities::Band(0.1 * Hz, 1.0 * Hz)); EXPECT_GT(tau_s(0), 0.0); EXPECT_GT(tau_s(1), 0.0); EXPECT_GT(tau_s(2), 0.0); @@ -99,7 +113,9 @@ TEST(Attenuation_ComputeTauSigma, DifferentFrequencyRanges) { // Wide range { - auto tau_s = compute_tau_sigma<3>(0.01, 1000.0); + auto tau_s = + compute_tau_sigma<3>(specfem::utilities::Band( + 0.01 * Hz, 1000.0 * Hz)); EXPECT_GT(tau_s(0), 0.0); EXPECT_GT(tau_s(1), 0.0); EXPECT_GT(tau_s(2), 0.0); @@ -111,8 +127,10 @@ TEST(Attenuation_ComputeTauSigma, DifferentFrequencyRanges) { TEST(Attenuation_ComputeTauSigma, TwoSLS) { constexpr type_real min_frequency = 0.1; constexpr type_real max_frequency = 100.0; + specfem::utilities::Band band(min_frequency * Hz, + max_frequency * Hz); - auto tau_s = compute_tau_sigma<2>(min_frequency, max_frequency); + auto tau_s = compute_tau_sigma<2>(band); EXPECT_EQ(tau_s.extent(0), 2); EXPECT_GT(tau_s(0), 0.0); diff --git a/tests/unit-tests/point/attenuation_tests.cpp b/tests/unit-tests/point/attenuation_tests.cpp index f8fea9f30..4e8ccecb9 100644 --- a/tests/unit-tests/point/attenuation_tests.cpp +++ b/tests/unit-tests/point/attenuation_tests.cpp @@ -122,14 +122,14 @@ TYPED_TEST(PointAttenuationTest, ValueConstructor2D) { rxz_val, rk_val); for (int i = 0; i < N; ++i) { + EXPECT_TRUE(specfem::utilities::is_close(att.kappa_relaxation_rate(i), + kappa_val(i))) + << "kappa_relaxation_rate(" << i + << "): " << ExpectedGot(kappa_val(i), att.kappa_relaxation_rate(i)); EXPECT_TRUE( - specfem::utilities::is_close(att.kappa_common_factor(i), kappa_val(i))) - << "kappa_common_factor(" << i - << "): " << ExpectedGot(kappa_val(i), att.kappa_common_factor(i)); - EXPECT_TRUE( - specfem::utilities::is_close(att.mu_common_factor(i), mu_val(i))) - << "mu_common_factor(" << i - << "): " << ExpectedGot(mu_val(i), att.mu_common_factor(i)); + specfem::utilities::is_close(att.mu_relaxation_rate(i), mu_val(i))) + << "mu_relaxation_rate(" << i + << "): " << ExpectedGot(mu_val(i), att.mu_relaxation_rate(i)); EXPECT_TRUE(specfem::utilities::is_close(att.Rxx(i), rxx_val(i))) << "Rxx(" << i << "): " << ExpectedGot(rxx_val(i), att.Rxx(i)); EXPECT_TRUE(specfem::utilities::is_close(att.Rxz(i), rxz_val(i))) @@ -169,10 +169,10 @@ TYPED_TEST(PointAttenuationTest, ValueConstructor3D) { rxz, ryz, rk); for (int i = 0; i < N; ++i) { + EXPECT_TRUE(specfem::utilities::is_close(att.kappa_relaxation_rate(i), + kappa_val(i))); EXPECT_TRUE( - specfem::utilities::is_close(att.kappa_common_factor(i), kappa_val(i))); - EXPECT_TRUE( - specfem::utilities::is_close(att.mu_common_factor(i), mu_val(i))); + specfem::utilities::is_close(att.mu_relaxation_rate(i), mu_val(i))); EXPECT_TRUE(specfem::utilities::is_close(att.Rxx(i), rxx(i))); EXPECT_TRUE(specfem::utilities::is_close(att.Ryy(i), ryy(i))); EXPECT_TRUE(specfem::utilities::is_close(att.Rxy(i), rxy(i))); @@ -522,8 +522,8 @@ TYPED_TEST(PointAttenuationTest, PrintMethod2D) { EXPECT_FALSE(s.empty()); EXPECT_NE(s.find("Attenuation Factors"), std::string::npos); - EXPECT_NE(s.find("kappa_common_factor"), std::string::npos); - EXPECT_NE(s.find("mu_common_factor"), std::string::npos); + EXPECT_NE(s.find("kappa_relaxation_rate"), std::string::npos); + EXPECT_NE(s.find("mu_relaxation_rate"), std::string::npos); EXPECT_NE(s.find("alpha_rk"), std::string::npos); EXPECT_NE(s.find("beta_rk"), std::string::npos); EXPECT_NE(s.find("gamma_rk"), std::string::npos);