From 240222e5ea5f74e5a2195d5405ee05f0b8fc8a5f Mon Sep 17 00:00:00 2001 From: Congyue Cui Date: Mon, 24 Feb 2025 17:14:18 -0500 Subject: [PATCH 1/4] refactor point_properties and property_container --- include/IO/property/reader.tpp | 18 +- include/IO/property/writer.tpp | 18 +- include/boundary_conditions/stacey/stacey.tpp | 20 +- .../acoustic/isotropic/frechet_derivative.hpp | 4 +- .../dim2/acoustic/isotropic/mass_matrix.tpp | 2 +- .../isotropic/properties_container.hpp | 218 +--------- .../medium/dim2/acoustic/isotropic/source.hpp | 2 +- .../medium/dim2/acoustic/isotropic/stress.hpp | 4 +- .../anisotropic/frechet_derivative.hpp | 14 +- .../anisotropic/properties_container.hpp | 341 +--------------- .../dim2/elastic/anisotropic/stress.hpp | 12 +- .../dim2/elastic/anisotropic/wavefield.hpp | 20 +- .../elastic/isotropic/frechet_derivative.hpp | 14 +- .../dim2/elastic/isotropic/mass_matrix.tpp | 6 +- .../isotropic/properties_container.hpp | 234 +---------- .../medium/dim2/elastic/isotropic/stress.hpp | 8 +- .../elastic/isotropic/stress_integrand.tpp | 10 +- .../dim2/elastic/isotropic/wavefield.hpp | 2 +- include/medium/properties_container.hpp | 159 ++++++++ include/point/properties.hpp | 382 +++++++----------- .../compute_wavefield/test_helper.hpp | 4 +- .../assembly/properties/properties.cpp | 245 ++++++----- 22 files changed, 563 insertions(+), 1174 deletions(-) diff --git a/include/IO/property/reader.tpp b/include/IO/property/reader.tpp index 60253ebb0..5adf72165 100644 --- a/include/IO/property/reader.tpp +++ b/include/IO/property/reader.tpp @@ -23,31 +23,19 @@ void specfem::IO::property_reader::read(specfem::compute::assembly { typename InputLibrary::Group elastic = file.openGroup("/ElasticIsotropic"); - elastic.openDataset("rho", properties.elastic_isotropic.h_rho).read(); - elastic.openDataset("mu", properties.elastic_isotropic.h_mu).read(); - elastic.openDataset("lambdaplus2mu", properties.elastic_isotropic.h_lambdaplus2mu).read(); + elastic.openDataset("data", properties.elastic_isotropic.h_data).read(); } { typename InputLibrary::Group elastic = file.openGroup("/ElasticAnisotropic"); - elastic.openDataset("rho", properties.elastic_anisotropic.h_rho).read(); - elastic.openDataset("c11", properties.elastic_anisotropic.h_c11).read(); - elastic.openDataset("c13", properties.elastic_anisotropic.h_c13).read(); - elastic.openDataset("c15", properties.elastic_anisotropic.h_c15).read(); - elastic.openDataset("c33", properties.elastic_anisotropic.h_c33).read(); - elastic.openDataset("c35", properties.elastic_anisotropic.h_c35).read(); - elastic.openDataset("c55", properties.elastic_anisotropic.h_c55).read(); - elastic.openDataset("c12", properties.elastic_anisotropic.h_c12).read(); - elastic.openDataset("c23", properties.elastic_anisotropic.h_c23).read(); - elastic.openDataset("c25", properties.elastic_anisotropic.h_c25).read(); + elastic.openDataset("data", properties.elastic_anisotropic.h_data).read(); } { typename InputLibrary::Group acoustic = file.openGroup("/Acoustic"); - acoustic.openDataset("rho_inverse", properties.acoustic_isotropic.h_rho_inverse).read(); - acoustic.openDataset("kappa", properties.acoustic_isotropic.h_kappa).read(); + acoustic.openDataset("data", properties.acoustic_isotropic.h_data).read(); } std::cout << "Properties read from " << input_folder << "/Properties" diff --git a/include/IO/property/writer.tpp b/include/IO/property/writer.tpp index 563a22428..b0682f0d5 100644 --- a/include/IO/property/writer.tpp +++ b/include/IO/property/writer.tpp @@ -57,9 +57,7 @@ void specfem::IO::property_writer::write(specfem::compute::assemb elastic.createDataset("X", x).write(); elastic.createDataset("Z", z).write(); - elastic.createDataset("rho", properties.elastic_isotropic.h_rho).write(); - elastic.createDataset("mu", properties.elastic_isotropic.h_mu).write(); - elastic.createDataset("lambdaplus2mu", properties.elastic_isotropic.h_lambdaplus2mu).write(); + elastic.createDataset("data", properties.elastic_isotropic.h_data).write(); } { @@ -86,16 +84,7 @@ void specfem::IO::property_writer::write(specfem::compute::assemb elastic.createDataset("X", x).write(); elastic.createDataset("Z", z).write(); - elastic.createDataset("rho", properties.elastic_anisotropic.h_rho).write(); - elastic.createDataset("c11", properties.elastic_anisotropic.h_c11).write(); - elastic.createDataset("c13", properties.elastic_anisotropic.h_c13).write(); - elastic.createDataset("c15", properties.elastic_anisotropic.h_c15).write(); - elastic.createDataset("c33", properties.elastic_anisotropic.h_c33).write(); - elastic.createDataset("c35", properties.elastic_anisotropic.h_c35).write(); - elastic.createDataset("c55", properties.elastic_anisotropic.h_c55).write(); - elastic.createDataset("c12", properties.elastic_anisotropic.h_c12).write(); - elastic.createDataset("c23", properties.elastic_anisotropic.h_c23).write(); - elastic.createDataset("c25", properties.elastic_anisotropic.h_c25).write(); + elastic.createDataset("data", properties.elastic_anisotropic.h_data).write(); } { @@ -120,8 +109,7 @@ void specfem::IO::property_writer::write(specfem::compute::assemb acoustic.createDataset("X", x).write(); acoustic.createDataset("Z", z).write(); - acoustic.createDataset("rho_inverse", properties.acoustic_isotropic.h_rho_inverse).write(); - acoustic.createDataset("kappa", properties.acoustic_isotropic.h_kappa).write(); + acoustic.createDataset("data", properties.acoustic_isotropic.h_data).write(); } assert(n_elastic_isotropic + n_elastic_anisotropic + n_acoustic == nspec); diff --git a/include/boundary_conditions/stacey/stacey.tpp b/include/boundary_conditions/stacey/stacey.tpp index 7f3010b0a..b84b4a2e5 100644 --- a/include/boundary_conditions/stacey/stacey.tpp +++ b/include/boundary_conditions/stacey/stacey.tpp @@ -57,7 +57,7 @@ impl_enforce_traction(const acoustic_type &, const isotropic_type &, // Apply Stacey boundary condition traction(0) += static_cast(-1.0) * factor * - property.rho_vpinverse * field.velocity(0); + property.rho_vpinverse() * field.velocity(0); // Do nothing return; @@ -101,7 +101,7 @@ impl_enforce_traction(const acoustic_type &, const isotropic_type &, // Apply Stacey boundary condition Kokkos::Experimental::where(mask, traction(0)) = traction(0) + static_cast(-1.0) * factor * - property.rho_vpinverse * field.velocity(0); + property.rho_vpinverse() * field.velocity(0); return; } @@ -146,8 +146,8 @@ impl_enforce_traction(const elastic_type &, const isotropic_type &, for (int icomp = 0; icomp < 2; ++icomp) { factor[icomp] = ((vn * dn(icomp) / (jacobian1d * jacobian1d)) * - (property.rho_vp - property.rho_vs)) + - field.velocity(icomp) * property.rho_vs; + (property.rho_vp() - property.rho_vs())) + + field.velocity(icomp) * property.rho_vs(); } traction(0) += static_cast(-1.0) * factor[0] * jacobian1d * @@ -203,8 +203,8 @@ impl_enforce_traction(const elastic_type &, const isotropic_type &, for (int icomp = 0; icomp < 2; ++icomp) { factor[icomp] = ((vn * dn(icomp) / (jacobian1d * jacobian1d)) * - (property.rho_vp - property.rho_vs)) + - field.velocity(icomp) * property.rho_vs; + (property.rho_vp() - property.rho_vs())) + + field.velocity(icomp) * property.rho_vs(); } Kokkos::Experimental::where(mask, traction(0)) = @@ -259,8 +259,8 @@ impl_enforce_traction(const elastic_type &, const anisotropic_type &, for (int icomp = 0; icomp < 2; ++icomp) { factor[icomp] = ((vn * dn(icomp) / (jacobian1d * jacobian1d)) * - (property.rho_vp - property.rho_vs)) + - field.velocity(icomp) * property.rho_vs; + (property.rho_vp() - property.rho_vs())) + + field.velocity(icomp) * property.rho_vs(); } traction(0) += static_cast(-1.0) * factor[0] * jacobian1d * @@ -316,8 +316,8 @@ impl_enforce_traction(const elastic_type &, const anisotropic_type &, for (int icomp = 0; icomp < 2; ++icomp) { factor[icomp] = ((vn * dn(icomp) / (jacobian1d * jacobian1d)) * - (property.rho_vp - property.rho_vs)) + - field.velocity(icomp) * property.rho_vs; + (property.rho_vp() - property.rho_vs())) + + field.velocity(icomp) * property.rho_vs(); } Kokkos::Experimental::where(mask, traction(0)) = diff --git a/include/medium/dim2/acoustic/isotropic/frechet_derivative.hpp b/include/medium/dim2/acoustic/isotropic/frechet_derivative.hpp index 25f48a344..59b2f2156 100644 --- a/include/medium/dim2/acoustic/isotropic/frechet_derivative.hpp +++ b/include/medium/dim2/acoustic/isotropic/frechet_derivative.hpp @@ -30,11 +30,11 @@ impl_compute_frechet_derivatives( const auto rho_kl = (adjoint_derivatives.du(0, 0) * backward_derivatives.du(0, 0) + adjoint_derivatives.du(1, 0) * backward_derivatives.du(1, 0)) * - properties.rho_inverse * dt; + properties.rho_inverse() * dt; const auto kappa_kl = specfem::algorithms::dot(adjoint_field.acceleration, backward_field.displacement) * - static_cast(1.0) / properties.kappa * dt; + static_cast(1.0) / properties.kappa() * dt; return { rho_kl, kappa_kl }; } diff --git a/include/medium/dim2/acoustic/isotropic/mass_matrix.tpp b/include/medium/dim2/acoustic/isotropic/mass_matrix.tpp index 4533bf91d..7af193841 100644 --- a/include/medium/dim2/acoustic/isotropic/mass_matrix.tpp +++ b/include/medium/dim2/acoustic/isotropic/mass_matrix.tpp @@ -14,5 +14,5 @@ specfem::medium::impl_mass_matrix_component( specfem::dimension::type::dim2, true, UseSIMD> &partial_derivatives) { return specfem::datatype::ScalarPointViewType( - partial_derivatives.jacobian / properties.kappa); + partial_derivatives.jacobian / properties.kappa()); } diff --git a/include/medium/dim2/acoustic/isotropic/properties_container.hpp b/include/medium/dim2/acoustic/isotropic/properties_container.hpp index 9c76c2b65..ee9381fec 100644 --- a/include/medium/dim2/acoustic/isotropic/properties_container.hpp +++ b/include/medium/dim2/acoustic/isotropic/properties_container.hpp @@ -9,211 +9,19 @@ namespace medium { template <> struct properties_container { - - constexpr static auto dimension = specfem::dimension::type::dim2; - constexpr static auto value_type = specfem::element::medium_tag::acoustic; - constexpr static auto property_type = - specfem::element::property_tag::isotropic; - - using ViewType = typename Kokkos::View; - - int nspec; ///< total number of acoustic spectral elements - int ngllz; ///< number of quadrature points in z dimension - int ngllx; ///< number of quadrature points in x dimension - ViewType rho_inverse; - ViewType::HostMirror h_rho_inverse; - ViewType kappa; - ViewType::HostMirror h_kappa; - - properties_container() = default; - - properties_container(const int nspec, const int ngllz, const int ngllx) - : nspec(nspec), ngllz(ngllz), ngllx(ngllx), - rho_inverse("specfem::compute::properties::rho_inverse", nspec, ngllz, - ngllx), - h_rho_inverse(Kokkos::create_mirror_view(rho_inverse)), - kappa("specfem::compute::properties::kappa", nspec, ngllz, ngllx), - h_kappa(Kokkos::create_mirror_view(kappa)) {} - - template < - typename PointProperties, - typename std::enable_if_t = 0> - KOKKOS_FORCEINLINE_FUNCTION void - load_device_properties(const specfem::point::index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - property.rho_inverse = rho_inverse(ispec, iz, ix); - property.kappa = kappa(ispec, iz, ix); - property.kappa_inverse = static_cast(1.0) / property.kappa; - property.rho_vpinverse = - sqrt(property.rho_inverse * property.kappa_inverse); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - KOKKOS_FORCEINLINE_FUNCTION void - load_device_properties(const specfem::point::simd_index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho_inverse) - .copy_from(&rho_inverse(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.kappa) - .copy_from(&kappa(ispec, iz, ix), tag_type()); - - property.kappa_inverse = static_cast(1.0) / property.kappa; - property.rho_vpinverse = - Kokkos::sqrt(property.rho_inverse * property.kappa_inverse); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void - load_host_properties(const specfem::point::index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - property.rho_inverse = h_rho_inverse(ispec, iz, ix); - property.kappa = h_kappa(ispec, iz, ix); - property.kappa_inverse = static_cast(1.0) / property.kappa; - property.rho_vpinverse = - sqrt(property.rho_inverse * property.kappa_inverse); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void - load_host_properties(const specfem::point::simd_index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho_inverse) - .copy_from(&h_rho_inverse(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.kappa) - .copy_from(&h_kappa(ispec, iz, ix), tag_type()); - - property.kappa_inverse = static_cast(1.0) / property.kappa; - property.rho_vpinverse = - Kokkos::sqrt(property.rho_inverse * property.kappa_inverse); - } - - void copy_to_device() { - Kokkos::deep_copy(rho_inverse, h_rho_inverse); - Kokkos::deep_copy(kappa, h_kappa); - } - - void copy_to_host() { - Kokkos::deep_copy(h_rho_inverse, rho_inverse); - Kokkos::deep_copy(h_kappa, kappa); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void assign(const specfem::point::index &index, - const PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - h_rho_inverse(ispec, iz, ix) = property.rho_inverse; - h_kappa(ispec, iz, ix) = property.kappa; - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void assign(const specfem::point::simd_index &index, - const PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho_inverse) - .copy_to(&h_rho_inverse(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.kappa) - .copy_to(&h_kappa(ispec, iz, ix), tag_type()); - } + specfem::element::property_tag::isotropic> + : public impl::impl_properties_container< + specfem::element::medium_tag::acoustic, + specfem::element::property_tag::isotropic, 2> { + using Base = + impl::impl_properties_container; + using Base::Base; + constexpr static int _counter = __COUNTER__; + + DEFINE_CONTAINER(rho_inverse) + DEFINE_CONTAINER(kappa) }; } // namespace medium diff --git a/include/medium/dim2/acoustic/isotropic/source.hpp b/include/medium/dim2/acoustic/isotropic/source.hpp index 6ba010adc..b5e7e43b3 100644 --- a/include/medium/dim2/acoustic/isotropic/source.hpp +++ b/include/medium/dim2/acoustic/isotropic/source.hpp @@ -27,7 +27,7 @@ KOKKOS_INLINE_FUNCTION auto impl_compute_source_contribution( result.acceleration(0) = point_source.stf(0) * point_source.lagrange_interpolant(0) / - point_properties.kappa; + point_properties.kappa(); return result; } diff --git a/include/medium/dim2/acoustic/isotropic/stress.hpp b/include/medium/dim2/acoustic/isotropic/stress.hpp index 895a43e52..38f2c5d6d 100644 --- a/include/medium/dim2/acoustic/isotropic/stress.hpp +++ b/include/medium/dim2/acoustic/isotropic/stress.hpp @@ -26,8 +26,8 @@ impl_compute_stress( specfem::datatype::VectorPointViewType T; - T(0, 0) = properties.rho_inverse * du(0, 0); - T(1, 0) = properties.rho_inverse * du(1, 0); + T(0, 0) = properties.rho_inverse() * du(0, 0); + T(1, 0) = properties.rho_inverse() * du(1, 0); return { T }; } diff --git a/include/medium/dim2/elastic/anisotropic/frechet_derivative.hpp b/include/medium/dim2/elastic/anisotropic/frechet_derivative.hpp index 2f83ec85c..97ee5a04b 100644 --- a/include/medium/dim2/elastic/anisotropic/frechet_derivative.hpp +++ b/include/medium/dim2/elastic/anisotropic/frechet_derivative.hpp @@ -79,13 +79,13 @@ impl_compute_frechet_derivatives( // Computing the rest of the integral. // rho from equation 14 - rho_kl = static_cast(-1.0) * properties.rho * dt * rho_kl; - c11_kl = static_cast(-1.0) * c11_kl * properties.c11 * dt; - c13_kl = static_cast(-1.0) * c13_kl * properties.c13 * dt; - c15_kl = static_cast(-1.0) * c15_kl * properties.c15 * dt; - c33_kl = static_cast(-1.0) * c33_kl * properties.c33 * dt; - c35_kl = static_cast(-1.0) * c35_kl * properties.c35 * dt; - c55_kl = static_cast(-1.0) * c55_kl * properties.c55 * dt; + rho_kl = static_cast(-1.0) * properties.rho() * dt * rho_kl; + c11_kl = static_cast(-1.0) * c11_kl * properties.c11() * dt; + c13_kl = static_cast(-1.0) * c13_kl * properties.c13() * dt; + c15_kl = static_cast(-1.0) * c15_kl * properties.c15() * dt; + c33_kl = static_cast(-1.0) * c33_kl * properties.c33() * dt; + c35_kl = static_cast(-1.0) * c35_kl * properties.c35() * dt; + c55_kl = static_cast(-1.0) * c55_kl * properties.c55() * dt; return { rho_kl, c11_kl, c13_kl, c15_kl, c33_kl, c35_kl, c55_kl }; diff --git a/include/medium/dim2/elastic/anisotropic/properties_container.hpp b/include/medium/dim2/elastic/anisotropic/properties_container.hpp index ca78ed052..b13079853 100644 --- a/include/medium/dim2/elastic/anisotropic/properties_container.hpp +++ b/include/medium/dim2/elastic/anisotropic/properties_container.hpp @@ -9,327 +9,26 @@ namespace medium { template <> struct properties_container { - - constexpr static auto dimension = specfem::dimension::type::dim2; - constexpr static auto value_type = specfem::element::medium_tag::elastic; - constexpr static auto property_type = - specfem::element::property_tag::anisotropic; - - using ViewType = typename Kokkos::View; - - int nspec; ///< total number of acoustic spectral elements - int ngllz; ///< number of quadrature points in z dimension - int ngllx; ///< number of quadrature points in x dimension - - ViewType rho; - ViewType::HostMirror h_rho; - ViewType c11; - ViewType::HostMirror h_c11; - ViewType c13; - ViewType::HostMirror h_c13; - ViewType c15; - ViewType::HostMirror h_c15; - ViewType c33; - ViewType::HostMirror h_c33; - ViewType c35; - ViewType::HostMirror h_c35; - ViewType c55; - ViewType::HostMirror h_c55; - ViewType c12; - ViewType::HostMirror h_c12; - ViewType c23; - ViewType::HostMirror h_c23; - ViewType c25; - ViewType::HostMirror h_c25; - - properties_container() = default; - - properties_container(const int nspec, const int ngllz, const int ngllx) - : nspec(nspec), ngllz(ngllz), ngllx(ngllx), - rho("specfem::compute::properties::rho", nspec, ngllz, ngllx), - h_rho(Kokkos::create_mirror_view(rho)), - c11("specfem::compute::properties::c11", nspec, ngllz, ngllx), - h_c11(Kokkos::create_mirror_view(c11)), - c12("specfem::compute::properties::c12", nspec, ngllz, ngllx), - h_c12(Kokkos::create_mirror_view(c12)), - c13("specfem::compute::properties::c13", nspec, ngllz, ngllx), - h_c13(Kokkos::create_mirror_view(c13)), - c15("specfem::compute::properties::c15", nspec, ngllz, ngllx), - h_c15(Kokkos::create_mirror_view(c15)), - c33("specfem::compute::properties::c33", nspec, ngllz, ngllx), - h_c33(Kokkos::create_mirror_view(c33)), - c35("specfem::compute::properties::c35", nspec, ngllz, ngllx), - h_c35(Kokkos::create_mirror_view(c35)), - c55("specfem::compute::properties::c55", nspec, ngllz, ngllx), - h_c55(Kokkos::create_mirror_view(c55)), - c23("specfem::compute::properties::c23", nspec, ngllz, ngllx), - h_c23(Kokkos::create_mirror_view(c23)), - c25("specfem::compute::properties::c25", nspec, ngllz, ngllx), - h_c25(Kokkos::create_mirror_view(c25)) {} - - template < - typename PointProperties, - typename std::enable_if_t = 0> - KOKKOS_FORCEINLINE_FUNCTION void - load_device_properties(const specfem::point::index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - property.rho = rho(ispec, iz, ix); - property.c11 = c11(ispec, iz, ix); - property.c12 = c12(ispec, iz, ix); - property.c13 = c13(ispec, iz, ix); - property.c15 = c15(ispec, iz, ix); - property.c33 = c33(ispec, iz, ix); - property.c35 = c35(ispec, iz, ix); - property.c55 = c55(ispec, iz, ix); - property.c23 = c23(ispec, iz, ix); - property.c25 = c25(ispec, iz, ix); - - property.rho_vp = sqrt(property.rho * property.c33); - property.rho_vs = sqrt(property.rho * property.c55); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - KOKKOS_FORCEINLINE_FUNCTION void - load_device_properties(const specfem::point::simd_index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho) - .copy_from(&rho(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c11) - .copy_from(&c11(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c12) - .copy_from(&c12(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c13) - .copy_from(&c13(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c15) - .copy_from(&c15(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c33) - .copy_from(&c33(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c35) - .copy_from(&c35(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c55) - .copy_from(&c55(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c23) - .copy_from(&c23(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c25) - .copy_from(&c25(ispec, iz, ix), tag_type()); - - property.rho_vp = Kokkos::sqrt(property.rho * property.c33); - property.rho_vs = Kokkos::sqrt(property.rho * property.c55); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void - load_host_properties(const specfem::point::index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - property.rho = h_rho(ispec, iz, ix); - property.c11 = h_c11(ispec, iz, ix); - property.c12 = h_c12(ispec, iz, ix); - property.c13 = h_c13(ispec, iz, ix); - property.c15 = h_c15(ispec, iz, ix); - property.c33 = h_c33(ispec, iz, ix); - property.c35 = h_c35(ispec, iz, ix); - property.c55 = h_c55(ispec, iz, ix); - property.c23 = h_c23(ispec, iz, ix); - property.c25 = h_c25(ispec, iz, ix); - - property.rho_vp = sqrt(property.rho * property.c33); - property.rho_vs = sqrt(property.rho * property.c55); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void - load_host_properties(const specfem::point::simd_index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho) - .copy_from(&h_rho(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c11) - .copy_from(&h_c11(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c13) - .copy_from(&h_c13(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c15) - .copy_from(&h_c15(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c33) - .copy_from(&h_c33(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c35) - .copy_from(&h_c35(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c55) - .copy_from(&h_c55(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c23) - .copy_from(&h_c23(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c25) - .copy_from(&h_c25(ispec, iz, ix), tag_type()); - - property.rho_vp = Kokkos::sqrt(property.rho * property.c33); - property.rho_vs = Kokkos::sqrt(property.rho * property.c55); - } - - void copy_to_device() { - Kokkos::deep_copy(rho, h_rho); - Kokkos::deep_copy(c11, h_c11); - Kokkos::deep_copy(c13, h_c13); - Kokkos::deep_copy(c15, h_c15); - Kokkos::deep_copy(c12, h_c12); - Kokkos::deep_copy(c33, h_c33); - Kokkos::deep_copy(c35, h_c35); - Kokkos::deep_copy(c55, h_c55); - Kokkos::deep_copy(c23, h_c23); - Kokkos::deep_copy(c25, h_c25); - } - - void copy_to_host() { - Kokkos::deep_copy(h_rho, rho); - Kokkos::deep_copy(h_c11, c11); - Kokkos::deep_copy(h_c13, c13); - Kokkos::deep_copy(h_c15, c15); - Kokkos::deep_copy(h_c12, c12); - Kokkos::deep_copy(h_c33, c33); - Kokkos::deep_copy(h_c35, c35); - Kokkos::deep_copy(h_c55, c55); - Kokkos::deep_copy(h_c23, c23); - Kokkos::deep_copy(h_c25, c25); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void assign(const specfem::point::index &index, - const PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - h_rho(ispec, iz, ix) = property.rho; - h_c11(ispec, iz, ix) = property.c11; - h_c13(ispec, iz, ix) = property.c13; - h_c15(ispec, iz, ix) = property.c15; - h_c12(ispec, iz, ix) = property.c12; - h_c33(ispec, iz, ix) = property.c33; - h_c35(ispec, iz, ix) = property.c35; - h_c55(ispec, iz, ix) = property.c55; - h_c23(ispec, iz, ix) = property.c23; - h_c25(ispec, iz, ix) = property.c25; - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void assign(const specfem::point::simd_index &index, - const PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho) - .copy_to(&h_rho(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c11) - .copy_to(&h_c11(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c13) - .copy_to(&h_c13(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c15) - .copy_to(&h_c15(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c12) - .copy_to(&h_c12(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c33) - .copy_to(&h_c33(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c35) - .copy_to(&h_c35(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c55) - .copy_to(&h_c55(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c23) - .copy_to(&h_c23(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.c25) - .copy_to(&h_c25(ispec, iz, ix), tag_type()); - } + specfem::element::property_tag::anisotropic> + : public impl::impl_properties_container< + specfem::element::medium_tag::elastic, + specfem::element::property_tag::anisotropic, 10> { + using Base = impl::impl_properties_container< + specfem::element::medium_tag::elastic, + specfem::element::property_tag::anisotropic, 10>; + using Base::Base; + constexpr static int _counter = __COUNTER__; + + DEFINE_CONTAINER(c11) + DEFINE_CONTAINER(c13) + DEFINE_CONTAINER(c15) + DEFINE_CONTAINER(c33) + DEFINE_CONTAINER(c35) + DEFINE_CONTAINER(c55) + DEFINE_CONTAINER(c12) + DEFINE_CONTAINER(c23) + DEFINE_CONTAINER(c25) + DEFINE_CONTAINER(rho) }; } // namespace medium diff --git a/include/medium/dim2/elastic/anisotropic/stress.hpp b/include/medium/dim2/elastic/anisotropic/stress.hpp index 39ba0382f..563cd2f45 100644 --- a/include/medium/dim2/elastic/anisotropic/stress.hpp +++ b/include/medium/dim2/elastic/anisotropic/stress.hpp @@ -28,16 +28,16 @@ impl_compute_stress( // P_SV case // sigma_xx - sigma_xx = properties.c11 * du(0, 0) + properties.c13 * du(1, 1) + - properties.c15 * (du(1, 0) + du(0, 1)); + sigma_xx = properties.c11() * du(0, 0) + properties.c13() * du(1, 1) + + properties.c15() * (du(1, 0) + du(0, 1)); // sigma_zz - sigma_zz = properties.c13 * du(0, 0) + properties.c33 * du(1, 1) + - properties.c35 * (du(1, 0) + du(0, 1)); + sigma_zz = properties.c13() * du(0, 0) + properties.c33() * du(1, 1) + + properties.c35() * (du(1, 0) + du(0, 1)); // sigma_xz - sigma_xz = properties.c15 * du(0, 0) + properties.c35 * du(1, 1) + - properties.c55 * (du(1, 0) + du(0, 1)); + sigma_xz = properties.c15() * du(0, 0) + properties.c35() * du(1, 1) + + properties.c55() * (du(1, 0) + du(0, 1)); specfem::datatype::VectorPointViewType T; diff --git a/include/medium/dim2/elastic/anisotropic/wavefield.hpp b/include/medium/dim2/elastic/anisotropic/wavefield.hpp index 4dfd2050b..b0f4ff6f9 100644 --- a/include/medium/dim2/elastic/anisotropic/wavefield.hpp +++ b/include/medium/dim2/elastic/anisotropic/wavefield.hpp @@ -66,7 +66,7 @@ KOKKOS_FUNCTION void impl_compute_wavefield( // cannot compute pressure for an anisotropic material if c12 or c23 // are zero - if (point_property.c12 < 1.e-7 || point_property.c23 < 1.e-7) { + if (point_property.c12() < 1.e-7 || point_property.c23() < 1.e-7) { Kokkos::abort("C_12 or C_23 are zero, cannot compute pressure. " "Check your material properties. Or, deactivate the " "pressure computation."); @@ -85,19 +85,19 @@ KOKKOS_FUNCTION void impl_compute_wavefield( // P_SV case // sigma_xx - const auto sigma_xx = point_property.c11 * du(0, 0) + - point_property.c13 * du(1, 1) + - point_property.c15 * (du(1, 0) + du(0, 1)); + const auto sigma_xx = point_property.c11() * du(0, 0) + + point_property.c13() * du(1, 1) + + point_property.c15() * (du(1, 0) + du(0, 1)); // sigma_zz - const auto sigma_zz = point_property.c13 * du(0, 0) + - point_property.c33 * du(1, 1) + - point_property.c35 * (du(1, 0) + du(0, 1)); + const auto sigma_zz = point_property.c13() * du(0, 0) + + point_property.c33() * du(1, 1) + + point_property.c35() * (du(1, 0) + du(0, 1)); // sigma_yy - const auto sigma_yy = point_property.c12 * du(0, 0) + - point_property.c23 * du(1, 1) + - point_property.c25 * (du(1, 0) + du(0, 1)); + const auto sigma_yy = point_property.c12() * du(0, 0) + + point_property.c23() * du(1, 1) + + point_property.c25() * (du(1, 0) + du(0, 1)); wavefield(iterator_index.ielement, index.iz, index.ix, 0) = -1.0 * (sigma_xx + sigma_zz + sigma_yy) / 3.0; diff --git a/include/medium/dim2/elastic/isotropic/frechet_derivative.hpp b/include/medium/dim2/elastic/isotropic/frechet_derivative.hpp index 3949495b1..8853fd435 100644 --- a/include/medium/dim2/elastic/isotropic/frechet_derivative.hpp +++ b/include/medium/dim2/elastic/isotropic/frechet_derivative.hpp @@ -34,7 +34,7 @@ impl_compute_frechet_derivatives( specfem::globals::simulation_wave == specfem::wave::sh, "Only P-SV and SH waves are supported."); - const auto kappa = properties.lambdaplus2mu - properties.mu; + const auto kappa = properties.lambdaplus2mu() - properties.mu(); if (specfem::globals::simulation_wave == specfem::wave::p_sv) { @@ -121,8 +121,8 @@ impl_compute_frechet_derivatives( // Finishing the kernels kappa_kl = static_cast(-1.0) * kappa * dt * kappa_kl; - mu_kl = static_cast(-2.0) * properties.mu * dt * mu_kl; - rho_kl = static_cast(-1.0) * properties.rho * dt * rho_kl; + mu_kl = static_cast(-2.0) * properties.mu() * dt * mu_kl; + rho_kl = static_cast(-1.0) * properties.rho() * dt * rho_kl; // rho' kernel, first term in Equation 20 const auto rhop_kl = rho_kl + kappa_kl + mu_kl; @@ -130,14 +130,14 @@ impl_compute_frechet_derivatives( // beta (shear wave) kernel, second term in Equation 20 const auto beta_kl = static_cast(2.0) * (mu_kl - static_cast(4.0 / 3.0) * - properties.mu / kappa * kappa_kl); + properties.mu() / kappa * kappa_kl); // alpha (compressional wave) kernel, third and last term in Eq. 20 // of Tromp et al 2005. const auto alpha_kl = static_cast(2.0) * (static_cast(1.0) + - static_cast(4.0 / 3.0) * properties.mu / kappa) * + static_cast(4.0 / 3.0) * properties.mu() / kappa) * kappa_kl; return { rho_kl, mu_kl, kappa_kl, rhop_kl, alpha_kl, beta_kl }; @@ -164,13 +164,13 @@ impl_compute_frechet_derivatives( */ const auto mu_kl = - static_cast(-2.0) * properties.mu * dt * + static_cast(-2.0) * properties.mu() * dt * static_cast(0.5) * // du#y_dx * duy_dx + (adjoint_derivatives.du(0, 0) * backward_derivatives.du(0, 0) + // du#y_dz * duy_dz adjoint_derivatives.du(1, 0) * backward_derivatives.du(1, 0)); - const auto rho_kl = static_cast(-1.0) * properties.rho * dt * + const auto rho_kl = static_cast(-1.0) * properties.rho() * dt * specfem::algorithms::dot(adjoint_field.acceleration, backward_field.displacement); const auto kappa_kl = decltype(mu_kl)(0.0); diff --git a/include/medium/dim2/elastic/isotropic/mass_matrix.tpp b/include/medium/dim2/elastic/isotropic/mass_matrix.tpp index f49f3d8f7..28b8ae7a0 100644 --- a/include/medium/dim2/elastic/isotropic/mass_matrix.tpp +++ b/include/medium/dim2/elastic/isotropic/mass_matrix.tpp @@ -15,11 +15,11 @@ specfem::medium::impl_mass_matrix_component( if constexpr (specfem::globals::simulation_wave == specfem::wave::p_sv) { return specfem::datatype::ScalarPointViewType( - partial_derivatives.jacobian * properties.rho, - partial_derivatives.jacobian * properties.rho); + partial_derivatives.jacobian * properties.rho(), + partial_derivatives.jacobian * properties.rho()); } else if constexpr (specfem::globals::simulation_wave == specfem::wave::sh) { return specfem::datatype::ScalarPointViewType( - partial_derivatives.jacobian * properties.rho, 0); + partial_derivatives.jacobian * properties.rho(), 0); } else { static_assert("Unknown wave type"); return specfem::datatype::ScalarPointViewType( diff --git a/include/medium/dim2/elastic/isotropic/properties_container.hpp b/include/medium/dim2/elastic/isotropic/properties_container.hpp index 4be8eda70..1a4e8bdcf 100644 --- a/include/medium/dim2/elastic/isotropic/properties_container.hpp +++ b/include/medium/dim2/elastic/isotropic/properties_container.hpp @@ -9,226 +9,20 @@ namespace medium { template <> struct properties_container { - - constexpr static auto dimension = specfem::dimension::type::dim2; - constexpr static auto value_type = specfem::element::medium_tag::elastic; - constexpr static auto property_type = - specfem::element::property_tag::isotropic; - - using ViewType = typename Kokkos::View; - - int nspec; ///< total number of acoustic spectral elements - int ngllz; ///< number of quadrature points in z dimension - int ngllx; ///< number of quadrature points in x dimension - ViewType rho; - ViewType::HostMirror h_rho; - ViewType mu; - ViewType::HostMirror h_mu; - ViewType lambdaplus2mu; - ViewType::HostMirror h_lambdaplus2mu; - - properties_container() = default; - - properties_container(const int nspec, const int ngllz, const int ngllx) - : nspec(nspec), ngllz(ngllz), ngllx(ngllx), - rho("specfem::compute::properties::rho", nspec, ngllz, ngllx), - h_rho(Kokkos::create_mirror_view(rho)), - mu("specfem::compute::properties::mu", nspec, ngllz, ngllx), - h_mu(Kokkos::create_mirror_view(mu)), - lambdaplus2mu("specfem::compute::properties::lambdaplus2mu", nspec, - ngllz, ngllx), - h_lambdaplus2mu(Kokkos::create_mirror_view(lambdaplus2mu)) {} - - template < - typename PointProperties, - typename std::enable_if_t = 0> - KOKKOS_FORCEINLINE_FUNCTION void - load_device_properties(const specfem::point::index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - property.rho = rho(ispec, iz, ix); - property.mu = mu(ispec, iz, ix); - property.lambdaplus2mu = lambdaplus2mu(ispec, iz, ix); - property.lambda = property.lambdaplus2mu - 2 * property.mu; - property.rho_vp = sqrt(property.rho * property.lambdaplus2mu); - property.rho_vs = sqrt(property.rho * property.mu); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - KOKKOS_FORCEINLINE_FUNCTION void - load_device_properties(const specfem::point::simd_index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho) - .copy_from(&rho(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.mu) - .copy_from(&mu(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.lambdaplus2mu) - .copy_from(&lambdaplus2mu(ispec, iz, ix), tag_type()); - - property.lambda = property.lambdaplus2mu - 2 * property.mu; - property.rho_vp = Kokkos::sqrt(property.rho * property.lambdaplus2mu); - property.rho_vs = Kokkos::sqrt(property.rho * property.mu); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void - load_host_properties(const specfem::point::index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - property.rho = h_rho(ispec, iz, ix); - property.mu = h_mu(ispec, iz, ix); - property.lambdaplus2mu = h_lambdaplus2mu(ispec, iz, ix); - property.lambda = property.lambdaplus2mu - 2 * property.mu; - property.rho_vp = sqrt(property.rho * property.lambdaplus2mu); - property.rho_vs = sqrt(property.rho * property.mu); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void - load_host_properties(const specfem::point::simd_index &index, - PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho) - .copy_from(&h_rho(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.mu) - .copy_from(&h_mu(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.lambdaplus2mu) - .copy_from(&h_lambdaplus2mu(ispec, iz, ix), tag_type()); - - property.lambda = property.lambdaplus2mu - 2 * property.mu; - property.rho_vp = Kokkos::sqrt(property.rho * property.lambdaplus2mu); - property.rho_vs = Kokkos::sqrt(property.rho * property.mu); - } - - void copy_to_device() { - Kokkos::deep_copy(rho, h_rho); - Kokkos::deep_copy(mu, h_mu); - Kokkos::deep_copy(lambdaplus2mu, h_lambdaplus2mu); - } - - void copy_to_host() { - Kokkos::deep_copy(h_rho, rho); - Kokkos::deep_copy(h_mu, mu); - Kokkos::deep_copy(h_lambdaplus2mu, lambdaplus2mu); - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void assign(const specfem::point::index &index, - const PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - h_rho(ispec, iz, ix) = property.rho; - h_mu(ispec, iz, ix) = property.mu; - h_lambdaplus2mu(ispec, iz, ix) = property.lambdaplus2mu; - } - - template < - typename PointProperties, - typename std::enable_if_t = 0> - inline void assign(const specfem::point::simd_index &index, - const PointProperties &property) const { - - static_assert(PointProperties::dimension == dimension, - "Dimension mismatch"); - static_assert(PointProperties::medium_tag == value_type, - "Medium tag mismatch"); - static_assert(PointProperties::property_tag == property_type, - "Property tag mismatch"); - - using simd = typename PointProperties::simd; - using mask_type = typename simd::mask_type; - using tag_type = typename simd::tag_type; - - const int ispec = index.ispec; - const int iz = index.iz; - const int ix = index.ix; - - mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); - - Kokkos::Experimental::where(mask, property.rho) - .copy_to(&h_rho(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.mu) - .copy_to(&h_mu(ispec, iz, ix), tag_type()); - Kokkos::Experimental::where(mask, property.lambdaplus2mu) - .copy_to(&h_lambdaplus2mu(ispec, iz, ix), tag_type()); - } + specfem::element::property_tag::isotropic> + : public impl::impl_properties_container< + specfem::element::medium_tag::elastic, + specfem::element::property_tag::isotropic, 3> { + using Base = + impl::impl_properties_container; + using Base::Base; + constexpr static int _counter = __COUNTER__; + + DEFINE_CONTAINER(lambdaplus2mu) + DEFINE_CONTAINER(mu) + DEFINE_CONTAINER(rho) }; } // namespace medium diff --git a/include/medium/dim2/elastic/isotropic/stress.hpp b/include/medium/dim2/elastic/isotropic/stress.hpp index a90f6cceb..e7088a7f9 100644 --- a/include/medium/dim2/elastic/isotropic/stress.hpp +++ b/include/medium/dim2/elastic/isotropic/stress.hpp @@ -29,13 +29,15 @@ impl_compute_stress( // P_SV case // sigma_xx - sigma_xx = properties.lambdaplus2mu * du(0, 0) + properties.lambda * du(1, 1); + sigma_xx = + properties.lambdaplus2mu() * du(0, 0) + properties.lambda() * du(1, 1); // sigma_zz - sigma_zz = properties.lambdaplus2mu * du(1, 1) + properties.lambda * du(0, 0); + sigma_zz = + properties.lambdaplus2mu() * du(1, 1) + properties.lambda() * du(0, 0); // sigma_xz - sigma_xz = properties.mu * (du(0, 1) + du(1, 0)); + sigma_xz = properties.mu() * (du(0, 1) + du(1, 0)); specfem::datatype::VectorPointViewType T; diff --git a/include/medium/dim2/elastic/isotropic/stress_integrand.tpp b/include/medium/dim2/elastic/isotropic/stress_integrand.tpp index e6d2f0d88..12f668a94 100644 --- a/include/medium/dim2/elastic/isotropic/stress_integrand.tpp +++ b/include/medium/dim2/elastic/isotropic/stress_integrand.tpp @@ -26,22 +26,22 @@ specfem::medium::impl_compute_stress_integrands( // P_SV case // sigma_xx sigma_xx = - properties.lambdaplus2mu * du(0, 0) + properties.lambda * du(1, 1); + properties.lambdaplus2mu() * du(0, 0) + properties.lambda() * du(1, 1); // sigma_zz sigma_zz = - properties.lambdaplus2mu * du(1, 1) + properties.lambda * du(0, 0); + properties.lambdaplus2mu() * du(1, 1) + properties.lambda() * du(0, 0); // sigma_xz - sigma_xz = properties.mu * (du(0, 1) + du(1, 0)); + sigma_xz = properties.mu() * (du(0, 1) + du(1, 0)); } else if (specfem::globals::simulation_wave == specfem::wave::sh) { // SH-case // sigma_xx - sigma_xx = properties.mu * du(0, 0); // would be sigma_xy in + sigma_xx = properties.mu() * du(0, 0); // would be sigma_xy in // CPU-version // sigma_xz - sigma_xz = properties.mu * du(1, 0); // sigma_zy + sigma_xz = properties.mu() * du(1, 0); // sigma_zy } specfem::datatype::VectorPointViewType F; diff --git a/include/medium/dim2/elastic/isotropic/wavefield.hpp b/include/medium/dim2/elastic/isotropic/wavefield.hpp index 5deb6b24e..7f4da8765 100644 --- a/include/medium/dim2/elastic/isotropic/wavefield.hpp +++ b/include/medium/dim2/elastic/isotropic/wavefield.hpp @@ -80,7 +80,7 @@ KOKKOS_FUNCTION void impl_compute_wavefield( // -1.0 * (sigma_xx + sigma_zz + sigma_yy) / 3.0; wavefield(iterator_index.ielement, index.iz, index.ix, 0) = -1.0 * - ((point_property.lambda + (2.0 / 3.0) * point_property.mu) * + ((point_property.lambda() + (2.0 / 3.0) * point_property.mu()) * (du(0, 0) + du(1, 1))); }); diff --git a/include/medium/properties_container.hpp b/include/medium/properties_container.hpp index 2b3c9d892..3fc831ce6 100644 --- a/include/medium/properties_container.hpp +++ b/include/medium/properties_container.hpp @@ -2,9 +2,166 @@ #include "enumerations/medium.hpp" +#define DEFINE_CONTAINER(prop) \ + constexpr static int i_##prop = __COUNTER__ - _counter - 1; \ + KOKKOS_INLINE_FUNCTION type_real &prop(const int &ispec, const int &iz, \ + const int &ix) const { \ + return Base::data(ispec, iz, ix, i_##prop); \ + } \ + KOKKOS_INLINE_FUNCTION type_real &h_##prop(const int &ispec, const int &iz, \ + const int &ix) const { \ + return Base::h_data(ispec, iz, ix, i_##prop); \ + } + namespace specfem { namespace medium { +namespace impl { +template +struct impl_properties_container { + using view_type = typename Kokkos::View; + constexpr static auto nprops = N; + constexpr static auto dimension = specfem::dimension::type::dim2; + constexpr static auto medium_tag = MediumTag; + constexpr static auto property_tag = PropertyTag; + + int nspec; ///< total number of acoustic spectral elements + int ngllz; ///< number of quadrature points in z dimension + int ngllx; ///< number of quadrature points in x dimension + + view_type data; + view_type::HostMirror h_data; + + impl_properties_container() = default; + + impl_properties_container(const int nspec, const int ngllz, const int ngllx) + : nspec(nspec), ngllz(ngllz), ngllx(ngllx), + data("specfem::benchmarks::properties::data", nspec, ngllz, ngllx, N), + h_data(Kokkos::create_mirror_view(data)) {} + +private: + template + KOKKOS_FORCEINLINE_FUNCTION void + load_properties(const specfem::point::index &index, + PointProperties &property, const ViewType &target) const { + + static_assert(PointProperties::dimension == dimension, + "Dimension mismatch"); + static_assert(PointProperties::medium_tag == medium_tag, + "Medium tag mismatch"); + static_assert(PointProperties::property_tag == property_tag, + "Property tag mismatch"); + + const int ispec = index.ispec; + const int iz = index.iz; + const int ix = index.ix; + + for (int i = 0; i < nprops; i++) { + property.data[i] = target(ispec, iz, ix, i); + } + + property.compute(); + } + template + KOKKOS_FORCEINLINE_FUNCTION void + load_properties(const specfem::point::simd_index &index, + PointProperties &property, const ViewType &target) const { + + static_assert(PointProperties::dimension == dimension, + "Dimension mismatch"); + static_assert(PointProperties::medium_tag == medium_tag, + "Medium tag mismatch"); + static_assert(PointProperties::property_tag == property_tag, + "Property tag mismatch"); + + using simd = typename PointProperties::simd; + using mask_type = typename simd::mask_type; + using tag_type = typename simd::tag_type; + + const int ispec = index.ispec; + const int iz = index.iz; + const int ix = index.ix; + + mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); + + for (int i = 0; i < nprops; i++) { + Kokkos::Experimental::where(mask, property.data[i]) + .copy_from(&target(ispec, iz, ix, i), tag_type()); + } + + property.compute(); + } + +public: + template + KOKKOS_FORCEINLINE_FUNCTION void + load_device_properties(const IndexType &index, + PointProperties &property) const { + load_properties(index, property, data); + } + + template + KOKKOS_FORCEINLINE_FUNCTION void + load_host_properties(const IndexType &index, + PointProperties &property) const { + load_properties(index, property, h_data); + } + + template + inline void assign(const specfem::point::index &index, + const PointProperties &property) const { + + static_assert(PointProperties::dimension == dimension, + "Dimension mismatch"); + static_assert(PointProperties::medium_tag == medium_tag, + "Medium tag mismatch"); + static_assert(PointProperties::property_tag == property_tag, + "Property tag mismatch"); + + const int ispec = index.ispec; + const int iz = index.iz; + const int ix = index.ix; + + for (int i = 0; i < nprops; i++) { + h_data(ispec, iz, ix, i) = property.data[i]; + } + } + + template + inline void assign(const specfem::point::simd_index &index, + const PointProperties &property) const { + + static_assert(PointProperties::dimension == dimension, + "Dimension mismatch"); + static_assert(PointProperties::medium_tag == medium_tag, + "Medium tag mismatch"); + static_assert(PointProperties::property_tag == property_tag, + "Property tag mismatch"); + + using simd = typename PointProperties::simd; + using mask_type = typename simd::mask_type; + using tag_type = typename simd::tag_type; + + const int ispec = index.ispec; + const int iz = index.iz; + const int ix = index.ix; + + mask_type mask([&, this](std::size_t lane) { return index.mask(lane); }); + + for (int i = 0; i < nprops; i++) { + Kokkos::Experimental::where(mask, property.data[i]) + .copy_to(&h_data(ispec, iz, ix, i), tag_type()); + } + } + + void copy_to_device() { Kokkos::deep_copy(data, h_data); } + + void copy_to_host() { Kokkos::deep_copy(h_data, data); } +}; +} // namespace impl + template struct properties_container { @@ -19,3 +176,5 @@ struct properties_container { #include "dim2/acoustic/isotropic/properties_container.hpp" #include "dim2/elastic/anisotropic/properties_container.hpp" #include "dim2/elastic/isotropic/properties_container.hpp" + +#undef DEFINE_CONTAINER diff --git a/include/point/properties.hpp b/include/point/properties.hpp index a64095232..fab5902dd 100644 --- a/include/point/properties.hpp +++ b/include/point/properties.hpp @@ -1,5 +1,4 @@ -#ifndef _POINT_PROPERTIES_HPP -#define _POINT_PROPERTIES_HPP +#pragma once #include "datatypes/simd.hpp" #include "enumerations/medium.hpp" @@ -7,6 +6,72 @@ namespace specfem { namespace point { +#define DEFINE_PROP(prop) \ + constexpr static int i_##prop = __COUNTER__ - _counter - 1; \ + KOKKOS_INLINE_FUNCTION value_type prop() const { \ + return Base::data[i_##prop]; \ + } \ + KOKKOS_INLINE_FUNCTION void prop(value_type val) { \ + Base::data[i_##prop] = val; \ + } + +namespace impl { +template struct impl_properties { + using simd = specfem::datatype::simd; ///< SIMD type + using value_type = typename simd::datatype; + constexpr static bool is_point_properties = true; + constexpr static auto nprops = N; + constexpr static auto nprops_extra = N_EX; + + value_type data[N_EX]; + + KOKKOS_FUNCTION + impl_properties() = default; + + /** + * @brief array constructor + * + */ + KOKKOS_FUNCTION + impl_properties(const value_type *value) { + for (int i = 0; i < N; ++i) { + data[i] = value[i]; + } + } + + /** + * @brief value constructor + * + */ + template = 0> + impl_properties(Args... args) : data{ args... } {} + + /** + * @brief Equality operator + * + */ + KOKKOS_FUNCTION + bool operator==(const impl_properties &rhs) const { + for (int i = 0; i < N_EX; ++i) { + if (data[i] != rhs.data[i]) { + return false; + } + } + return true; + } + + /** + * @brief Inequality operator + * + */ + KOKKOS_FUNCTION + bool operator!=(const impl_properties &rhs) const { + return !(*this == rhs); + } +}; +} // namespace impl + /** * @brief Store properties of the medium at a quadrature point * @@ -28,216 +93,110 @@ struct properties; template struct properties { + specfem::element::property_tag::isotropic, UseSIMD> + : public impl::impl_properties<3, 6, UseSIMD> { /** * @name Typedefs * */ ///@{ - constexpr static bool is_point_properties = true; - using simd = specfem::datatype::simd; ///< SIMD type + using Base = impl::impl_properties<3, 6, UseSIMD>; constexpr static auto dimension = specfem::dimension::type::dim2; constexpr static auto medium_tag = specfem::element::medium_tag::elastic; constexpr static auto property_tag = specfem::element::property_tag::isotropic; - using value_type = - typename simd::datatype; ///< Value type to store properties - ///@} - value_type mu; ///< shear modulus @f$ \mu @f$ - value_type rho; ///< density @f$ \rho @f$ - - value_type rho_vp; ///< P-wave velocity @f$ \rho v_p @f$ - value_type rho_vs; ///< S-wave velocity @f$ \rho v_s @f$ - value_type lambda; ///< Lame's parameter @f$ \lambda @f$ - value_type lambdaplus2mu; ///< Lame's parameter @f$ \lambda + 2\mu @f$ - -private: - KOKKOS_FUNCTION - properties(const value_type &lambdaplus2mu, const value_type &mu, - const value_type &rho, std::false_type) - : lambdaplus2mu(lambdaplus2mu), mu(mu), rho(rho), - rho_vp(sqrt(rho * lambdaplus2mu)), rho_vs(sqrt(rho * mu)), - lambda(lambdaplus2mu - 2.0 * mu) {} - - KOKKOS_FUNCTION - properties(const value_type &lambdaplus2mu, const value_type &mu, - const value_type &rho, std::true_type) - : lambdaplus2mu(lambdaplus2mu), mu(mu), rho(rho), - rho_vp(Kokkos::sqrt(rho * lambdaplus2mu)), - rho_vs(Kokkos::sqrt(rho * mu)), - lambda(lambdaplus2mu - (static_cast(2.0)) * mu) {} - -public: - /** - * @name Constructors - * - */ - ///@{ - /** - * @brief Default constructor - * - */ + using simd = specfem::datatype::simd; + using value_type = typename simd::datatype; + constexpr static int _counter = __COUNTER__; + ///@} KOKKOS_FUNCTION properties() = default; - /** - * @brief Construct a new properties object - * - * @param lambdaplus2mu @f$ \lambda + 2\mu @f$ - * @param mu @f$ \mu @f$ - * @param rho @f$ \rho @f$ - */ - KOKKOS_FUNCTION - properties(const value_type &lambdaplus2mu, const value_type &mu, - const value_type &rho) - : properties(lambdaplus2mu, mu, rho, - std::integral_constant{}) {} - - /** - * @brief single value constructor - * - */ KOKKOS_FUNCTION - properties(const value_type value) - : properties(value, value, value, - std::integral_constant{}) {} - ///@} + properties(const value_type *value) : Base(value) { compute(); } - /** - * @brief Equality operator - * - */ - KOKKOS_FUNCTION - bool operator==(const properties &rhs) const { - return rho == rhs.rho && mu == rhs.mu && lambdaplus2mu == rhs.lambdaplus2mu; + template = 0> + KOKKOS_FUNCTION properties(Args... args) : Base(args...) { + compute(); } - /** - * @brief Inequality operator - * - */ - KOKKOS_FUNCTION - bool operator!=(const properties &rhs) const { return !(*this == rhs); } + DEFINE_PROP(lambdaplus2mu) ///< Lame's parameter @f$ \lambda + 2\mu @f$ + DEFINE_PROP(mu) ///< shear modulus @f$ \mu @f$ + DEFINE_PROP(rho) ///< density @f$ \rho @f$ - KOKKOS_FUNCTION - bool operator==(const value_type value) { - return rho == value && mu == value && lambdaplus2mu == value; - } + DEFINE_PROP(rho_vp) ///< P-wave velocity @f$ \rho v_p @f$ + DEFINE_PROP(rho_vs) ///< S-wave velocity @f$ \rho v_s @f$ + DEFINE_PROP(lambda) ///< Lame's parameter @f$ \lambda @f$ - KOKKOS_FUNCTION - bool operator!=(const value_type value) { return !(*this == value); } + KOKKOS_INLINE_FUNCTION + void compute() { + rho_vp(Kokkos::sqrt(rho() * lambdaplus2mu())); + rho_vs(Kokkos::sqrt(rho() * mu())); + lambda(lambdaplus2mu() - (static_cast(2.0)) * mu()); + } }; template struct properties { + specfem::element::property_tag::anisotropic, UseSIMD> + : public impl::impl_properties<10, 12, UseSIMD> { /** * @name Typedefs * */ ///@{ - constexpr static bool is_point_properties = true; - using simd = specfem::datatype::simd; ///< SIMD type + using Base = impl::impl_properties<10, 12, UseSIMD>; constexpr static auto dimension = specfem::dimension::type::dim2; constexpr static auto medium_tag = specfem::element::medium_tag::elastic; constexpr static auto property_tag = specfem::element::property_tag::anisotropic; - using value_type = - typename simd::datatype; ///< Value type to store properties - ///@} - /** - * @name Elastic constants - * - */ - ///@{ - value_type c11; ///< @f$ c_{11} @f$ - value_type c13; ///< @f$ c_{13} @f$ - value_type c15; ///< @f$ c_{15} @f$ - value_type c33; ///< @f$ c_{33} @f$ - value_type c35; ///< @f$ c_{35} @f$ - value_type c55; ///< @f$ c_{55} @f$ - value_type c12; ///< @f$ c_{12} @f$ - value_type c23; ///< @f$ c_{23} @f$ - value_type c25; ///< @f$ c_{25} @f$ + using simd = specfem::datatype::simd; + using value_type = typename simd::datatype; + constexpr static int _counter = __COUNTER__; ///@} - - value_type rho; ///< Density @f$ \rho @f$ - value_type rho_vp; ///< P-wave velocity @f$ \rho v_p @f$ - value_type rho_vs; ///< S-wave velocity @f$ \rho v_s @f$ - -private: - KOKKOS_FUNCTION - properties(const value_type &c11, const value_type &c13, - const value_type &c15, const value_type &c33, - const value_type &c35, const value_type &c55, - const value_type &c12, const value_type &c23, - const value_type &c25, const value_type &rho, std::false_type) - : c11(c11), c13(c13), c15(c15), c33(c33), c35(c35), c55(c55), c12(c12), - c23(c23), c25(c25), rho(rho), rho_vp(sqrt(rho * c33)), - rho_vs(sqrt(rho * c55)) {} - - KOKKOS_FUNCTION - properties(const value_type &c11, const value_type &c13, - const value_type &c15, const value_type &c33, - const value_type &c35, const value_type &c55, - const value_type &c12, const value_type &c23, - const value_type &c25, const value_type &rho, std::true_type) - : c11(c11), c13(c13), c15(c15), c33(c33), c35(c35), c55(c55), c12(c12), - c23(c23), c25(c25), rho(rho), rho_vp(Kokkos::sqrt(rho * c33)), - rho_vs(Kokkos::sqrt(rho * c55)) {} - -public: KOKKOS_FUNCTION properties() = default; KOKKOS_FUNCTION - properties(const value_type &c11, const value_type &c13, - const value_type &c15, const value_type &c33, - const value_type &c35, const value_type &c55, - const value_type &c12, const value_type &c23, - const value_type &c25, const type_real &rho) - : properties(c11, c13, c15, c33, c35, c55, c12, c23, c25, rho, - std::integral_constant{}) {} - - /** - * @brief single value constructor - * - */ - KOKKOS_FUNCTION - properties(const value_type value) - : properties(value, value, value, value, value, value, value, value, - value, value, std::integral_constant{}) {} + properties(const value_type *value) : Base(value) { compute(); } - /** - * @brief Equality operator - * - */ - KOKKOS_FUNCTION - bool operator==(const properties &rhs) const { - return rho == rhs.rho && c11 == rhs.c11 && c13 == rhs.c13 && - c15 == rhs.c15 && c33 == rhs.c33 && c35 == rhs.c35 && c55 == rhs.c55; + template = 0> + KOKKOS_FUNCTION properties(Args... args) : Base(args...) { + compute(); } /** - * @brief Inequality operator + * @name Elastic constants * */ - KOKKOS_FUNCTION - bool operator!=(const properties &rhs) const { return !(*this == rhs); } + ///@{ + DEFINE_PROP(c11) ///< @f$ c_{11} @f$ + DEFINE_PROP(c13) ///< @f$ c_{13} @f$ + DEFINE_PROP(c15) ///< @f$ c_{15} @f$ + DEFINE_PROP(c33) ///< @f$ c_{33} @f$ + DEFINE_PROP(c35) ///< @f$ c_{35} @f$ + DEFINE_PROP(c55) ///< @f$ c_{55} @f$ + DEFINE_PROP(c12) ///< @f$ c_{12} @f$ + DEFINE_PROP(c23) ///< @f$ c_{23} @f$ + DEFINE_PROP(c25) ///< @f$ c_{25} @f$ + ///@} - KOKKOS_FUNCTION - bool operator==(const value_type value) { - return rho == value && c11 == value && c13 == value && c15 == value && - c33 == value && c35 == value && c55 == value; - } + DEFINE_PROP(rho) ///< Density @f$ \rho @f$ + DEFINE_PROP(rho_vp) ///< P-wave velocity @f$ \rho v_p @f$ + DEFINE_PROP(rho_vs) ///< S-wave velocity @f$ \rho v_s @f$ - KOKKOS_FUNCTION - bool operator!=(const value_type value) { return !(*this == value); } + KOKKOS_INLINE_FUNCTION + void compute() { + rho_vp(Kokkos::sqrt(rho() * c33())); + rho_vs(Kokkos::sqrt(rho() * c55())); + } }; /** @@ -248,103 +207,50 @@ struct properties struct properties { + specfem::element::property_tag::isotropic, UseSIMD> + : public impl::impl_properties<2, 4, UseSIMD> { /** * @name Typedefs * */ ///@{ - constexpr static bool is_point_properties = true; - using simd = specfem::datatype::simd; ///< SIMD type + using Base = impl::impl_properties<2, 4, UseSIMD>; constexpr static auto dimension = specfem::dimension::type::dim2; constexpr static auto medium_tag = specfem::element::medium_tag::acoustic; constexpr static auto property_tag = specfem::element::property_tag::isotropic; - using value_type = - typename simd::datatype; ///< Value type to store properties + using simd = specfem::datatype::simd; + using value_type = typename simd::datatype; + constexpr static int _counter = __COUNTER__; ///@} - - value_type kappa_inverse; ///< @f$ \frac{1}{\lambda + 2\mu} @f$ - value_type rho_inverse; ///< @f$ \frac{1}{\rho} @f$ - value_type kappa; ///< Bulk modulus @f$ \kappa @f$ - - value_type rho_vpinverse; ///< @f$ \frac{1}{\rho v_p} @f$ - -private: - KOKKOS_FUNCTION - properties(const value_type &rho_inverse, const value_type &kappa, - std::false_type) - : kappa_inverse(1.0 / kappa), rho_inverse(rho_inverse), kappa(kappa), - rho_vpinverse(sqrt(rho_inverse * kappa_inverse)) {} - - KOKKOS_FUNCTION - properties(const value_type &rho_inverse, const value_type &kappa, - std::true_type) - : kappa_inverse((static_cast(1.0)) / kappa), - rho_inverse(rho_inverse), kappa(kappa), - rho_vpinverse(Kokkos::sqrt(rho_inverse * kappa_inverse)) {} - -public: - /** - * @name Constructors - * - */ - ///@{ - - /** - * @brief Default constructor - * - */ KOKKOS_FUNCTION properties() = default; - /** - * @brief Construct a new properties object - * - * @param rho_inverse @f$ \frac{1}{\rho} @f$ - * @param kappa Bulk modulus @f$ \kappa @f$ - */ - KOKKOS_FUNCTION - properties(const value_type &rho_inverse, const value_type &kappa) - : properties(rho_inverse, kappa, - std::integral_constant{}) {} - /** - * @brief single value constructor - * - */ KOKKOS_FUNCTION - properties(const value_type value) - : properties(value, value, std::integral_constant{}) {} - ///@} + properties(const value_type *value) : Base(value) { compute(); } - /** - * @brief Equality operator - * - */ - KOKKOS_FUNCTION - bool operator==(const properties &rhs) const { - return rho_inverse == rhs.rho_inverse && kappa == rhs.kappa; + template = 0> + KOKKOS_FUNCTION properties(Args... args) : Base(args...) { + compute(); } - /** - * @brief Inequality operator - * - */ - KOKKOS_FUNCTION - bool operator!=(const properties &rhs) const { return !(*this == rhs); } + DEFINE_PROP(rho_inverse) ///< @f$ \frac{1}{\rho} @f$ + DEFINE_PROP(kappa) ///< Bulk modulus @f$ \kappa @f$ + DEFINE_PROP(kappa_inverse) ///< @f$ \frac{1}{\lambda + 2\mu} @f$ + DEFINE_PROP(rho_vpinverse) ///< @f$ \frac{1}{\rho v_p} @f$ - KOKKOS_FUNCTION - bool operator==(const value_type value) { - return rho_inverse == value && kappa == value; + KOKKOS_INLINE_FUNCTION + void compute() { + kappa_inverse((static_cast(1.0)) / kappa()); + rho_vpinverse(Kokkos::sqrt(rho_inverse() * kappa_inverse())); } - - KOKKOS_FUNCTION - bool operator!=(const value_type value) { return !(*this == value); } }; +#undef DEFINE_PROP_CONSTRUCTORS +#undef DEFINE_PROP + } // namespace point } // namespace specfem - -#endif diff --git a/tests/unit-tests/assembly/compute_wavefield/test_helper.hpp b/tests/unit-tests/assembly/compute_wavefield/test_helper.hpp index 52d45d06f..58ad19e6e 100644 --- a/tests/unit-tests/assembly/compute_wavefield/test_helper.hpp +++ b/tests/unit-tests/assembly/compute_wavefield/test_helper.hpp @@ -111,7 +111,7 @@ class test_helper 1.0e-4) { diff --git a/tests/unit-tests/assembly/properties/properties.cpp b/tests/unit-tests/assembly/properties/properties.cpp index c047bb029..abe828ac9 100644 --- a/tests/unit-tests/assembly/properties/properties.cpp +++ b/tests/unit-tests/assembly/properties/properties.cpp @@ -38,9 +38,9 @@ std::string get_error_message( std::ostringstream message; error_message_header(message, value, mode); - message << "\t\trho = " << point_property.rho << "\n"; - message << "\t\tmu = " << point_property.mu << "\n"; - message << "\t\tkappa = " << point_property.lambdaplus2mu << "\n"; + message << "\t\trho = " << point_property.rho() << "\n"; + message << "\t\tmu = " << point_property.mu() << "\n"; + message << "\t\tlambdaplus2mu = " << point_property.lambdaplus2mu() << "\n"; return message.str(); } @@ -54,13 +54,13 @@ std::string get_error_message( std::ostringstream message; error_message_header(message, value, mode); - message << "\t\trho = " << point_property.rho << "\n"; - message << "\t\tc11 = " << point_property.c11 << "\n"; - message << "\t\tc13 = " << point_property.c13 << "\n"; - message << "\t\tc15 = " << point_property.c15 << "\n"; - message << "\t\tc33 = " << point_property.c33 << "\n"; - message << "\t\tc35 = " << point_property.c35 << "\n"; - message << "\t\tc55 = " << point_property.c55 << "\n"; + message << "\t\trho = " << point_property.rho() << "\n"; + message << "\t\tc11 = " << point_property.c11() << "\n"; + message << "\t\tc13 = " << point_property.c13() << "\n"; + message << "\t\tc15 = " << point_property.c15() << "\n"; + message << "\t\tc33 = " << point_property.c33() << "\n"; + message << "\t\tc35 = " << point_property.c35() << "\n"; + message << "\t\tc55 = " << point_property.c55() << "\n"; return message.str(); } @@ -74,8 +74,8 @@ std::string get_error_message( std::ostringstream message; error_message_header(message, value, mode); - message << "\t\trho_inverse = " << point_property.rho_inverse << "\n"; - message << "\t\tkappa = " << point_property.kappa << "\n"; + message << "\t\trho_inverse = " << point_property.rho_inverse() << "\n"; + message << "\t\tkappa = " << point_property.kappa() << "\n"; return message.str(); } @@ -112,10 +112,10 @@ get_point_property(const int ispec, const int iz, const int ix, specfem::element::property_tag::isotropic, false> point_property; - point_property.rho = elastic_isotropic.h_rho(ispec_l, iz, ix); - point_property.mu = elastic_isotropic.h_mu(ispec_l, iz, ix); - point_property.lambdaplus2mu = - elastic_isotropic.h_lambdaplus2mu(ispec_l, iz, ix); + point_property.rho(elastic_isotropic.h_rho(ispec_l, iz, ix)); + point_property.mu(elastic_isotropic.h_mu(ispec_l, iz, ix)); + point_property.lambdaplus2mu( + elastic_isotropic.h_lambdaplus2mu(ispec_l, iz, ix)); return point_property; } @@ -134,9 +134,9 @@ get_point_property( specfem::element::property_tag::isotropic, false> point_property_l; - point_property_l.rho = point_property.rho[lane]; - point_property_l.mu = point_property.mu[lane]; - point_property_l.lambdaplus2mu = point_property.lambdaplus2mu[lane]; + point_property_l.rho(point_property.rho()[lane]); + point_property_l.mu(point_property.mu()[lane]); + point_property_l.lambdaplus2mu(point_property.lambdaplus2mu()[lane]); return point_property_l; } @@ -157,13 +157,16 @@ get_point_property(const int ispec, const int iz, const int ix, specfem::element::property_tag::anisotropic, false> point_property; - point_property.rho = elastic_anisotropic.h_rho(ispec_l, iz, ix); - point_property.c11 = elastic_anisotropic.h_c11(ispec_l, iz, ix); - point_property.c13 = elastic_anisotropic.h_c13(ispec_l, iz, ix); - point_property.c15 = elastic_anisotropic.h_c15(ispec_l, iz, ix); - point_property.c33 = elastic_anisotropic.h_c33(ispec_l, iz, ix); - point_property.c35 = elastic_anisotropic.h_c35(ispec_l, iz, ix); - point_property.c55 = elastic_anisotropic.h_c55(ispec_l, iz, ix); + point_property.rho(elastic_anisotropic.h_rho(ispec_l, iz, ix)); + point_property.c11(elastic_anisotropic.h_c11(ispec_l, iz, ix)); + point_property.c13(elastic_anisotropic.h_c13(ispec_l, iz, ix)); + point_property.c15(elastic_anisotropic.h_c15(ispec_l, iz, ix)); + point_property.c33(elastic_anisotropic.h_c33(ispec_l, iz, ix)); + point_property.c35(elastic_anisotropic.h_c35(ispec_l, iz, ix)); + point_property.c55(elastic_anisotropic.h_c55(ispec_l, iz, ix)); + point_property.c12(elastic_anisotropic.h_c12(ispec_l, iz, ix)); + point_property.c23(elastic_anisotropic.h_c23(ispec_l, iz, ix)); + point_property.c25(elastic_anisotropic.h_c25(ispec_l, iz, ix)); return point_property; } @@ -182,13 +185,16 @@ get_point_property( specfem::element::property_tag::anisotropic, false> point_property_l; - point_property_l.rho = point_property.rho[lane]; - point_property_l.c11 = point_property.c11[lane]; - point_property_l.c13 = point_property.c13[lane]; - point_property_l.c15 = point_property.c15[lane]; - point_property_l.c33 = point_property.c33[lane]; - point_property_l.c35 = point_property.c35[lane]; - point_property_l.c55 = point_property.c55[lane]; + point_property_l.rho(point_property.rho()[lane]); + point_property_l.c11(point_property.c11()[lane]); + point_property_l.c13(point_property.c13()[lane]); + point_property_l.c15(point_property.c15()[lane]); + point_property_l.c33(point_property.c33()[lane]); + point_property_l.c35(point_property.c35()[lane]); + point_property_l.c55(point_property.c55()[lane]); + point_property_l.c12(point_property.c12()[lane]); + point_property_l.c23(point_property.c23()[lane]); + point_property_l.c25(point_property.c25()[lane]); return point_property_l; } @@ -209,9 +215,8 @@ get_point_property(const int ispec, const int iz, const int ix, specfem::element::property_tag::isotropic, false> point_property; - point_property.rho_inverse = - acoustic_isotropic.h_rho_inverse(ispec_l, iz, ix); - point_property.kappa = acoustic_isotropic.h_kappa(ispec_l, iz, ix); + point_property.rho_inverse(acoustic_isotropic.h_rho_inverse(ispec_l, iz, ix)); + point_property.kappa(acoustic_isotropic.h_kappa(ispec_l, iz, ix)); return point_property; } @@ -230,8 +235,8 @@ get_point_property( specfem::element::property_tag::isotropic, false> point_property_l; - point_property_l.rho_inverse = point_property.rho_inverse[lane]; - point_property_l.kappa = point_property.kappa[lane]; + point_property_l.rho_inverse(point_property.rho_inverse()[lane]); + point_property_l.kappa(point_property.kappa()[lane]); return point_property_l; } @@ -240,13 +245,13 @@ template void check_eq( const typename specfem::datatype::simd::datatype &p1, const typename specfem::datatype::simd::datatype &p2, - const int &n_simd_elements) { + const int &n_simd_elements, const std::string &msg) { if constexpr (using_simd) { for (int i = 0; i < n_simd_elements; i++) { if (p1[i] != p2[i]) { std::ostringstream message; - message << "\n \t Error in function load_on_host"; + message << "\n \t Error in function load_on_host (SIMD) | " << msg; message << "\n\t Expected: " << p1[i]; message << "\n\t Got: " << p2[i]; @@ -257,7 +262,7 @@ void check_eq( if (p1 != p2) { std::ostringstream message; - message << "\n \t Error in function load_on_host"; + message << "\n \t Error in function load_on_host | " << msg; message << "\n\t Expected: " << p1; message << "\n\t Got: " << p2; @@ -284,19 +289,20 @@ void check_point_properties( specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, specfem::element::property_tag::isotropic, using_simd> &p2, const int &n_simd_elements) { - check_eq(p1.rho, p2.rho, n_simd_elements); - check_eq(p1.mu, p2.mu, n_simd_elements); - check_eq(p1.lambdaplus2mu, p2.lambdaplus2mu, n_simd_elements); - check_eq(p1.lambda, - p2.lambdaplus2mu - + check_eq(p1.rho(), p2.rho(), n_simd_elements, "rho"); + check_eq(p1.mu(), p2.mu(), n_simd_elements, ".mu"); + check_eq(p1.lambdaplus2mu(), p2.lambdaplus2mu(), n_simd_elements, + "lambdaplus2mu"); + check_eq(p1.lambda(), + p2.lambdaplus2mu() - (static_cast::datatype>(2.0)) * - p2.mu, - n_simd_elements); - check_eq(p1.rho_vp, Kokkos::sqrt(p2.rho * p2.lambdaplus2mu), - n_simd_elements); - check_eq(p1.rho_vs, Kokkos::sqrt(p2.rho * p2.mu), - n_simd_elements); + p2.mu(), + n_simd_elements, "lambda"); + check_eq(p1.rho_vp(), Kokkos::sqrt(p2.rho() * p2.lambdaplus2mu()), + n_simd_elements, "rho_vp"); + check_eq(p1.rho_vs(), Kokkos::sqrt(p2.rho() * p2.mu()), + n_simd_elements, "rho_vp"); } template @@ -308,17 +314,17 @@ void check_point_properties( specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, specfem::element::property_tag::anisotropic, using_simd> &p2, const int &n_simd_elements) { - check_eq(p1.rho, p2.rho, n_simd_elements); - check_eq(p1.c11, p2.c11, n_simd_elements); - check_eq(p1.c13, p2.c13, n_simd_elements); - check_eq(p1.c15, p2.c15, n_simd_elements); - check_eq(p1.c33, p2.c33, n_simd_elements); - check_eq(p1.c35, p2.c35, n_simd_elements); - check_eq(p1.c55, p2.c55, n_simd_elements); - check_eq(p1.rho_vp, Kokkos::sqrt(p2.rho * p2.c33), - n_simd_elements); - check_eq(p1.rho_vs, Kokkos::sqrt(p2.rho * p2.c55), - n_simd_elements); + check_eq(p1.rho(), p2.rho(), n_simd_elements, "rho"); + check_eq(p1.c11(), p2.c11(), n_simd_elements, "c11"); + check_eq(p1.c13(), p2.c13(), n_simd_elements, "c13"); + check_eq(p1.c15(), p2.c15(), n_simd_elements, "c15"); + check_eq(p1.c33(), p2.c33(), n_simd_elements, "c33"); + check_eq(p1.c35(), p2.c35(), n_simd_elements, "c35"); + check_eq(p1.c55(), p2.c55(), n_simd_elements, "c55"); + check_eq(p1.rho_vp(), Kokkos::sqrt(p2.rho() * p2.c33()), + n_simd_elements, "rho_vp"); + check_eq(p1.rho_vs(), Kokkos::sqrt(p2.rho() * p2.c55()), + n_simd_elements, "rho_vs"); } template @@ -330,18 +336,23 @@ void check_point_properties( specfem::dimension::type::dim2, specfem::element::medium_tag::acoustic, specfem::element::property_tag::isotropic, using_simd> &p2, const int &n_simd_elements) { - check_eq(p1.rho_inverse, p2.rho_inverse, n_simd_elements); - check_eq(p1.kappa, p2.kappa, n_simd_elements); + check_eq(p1.rho_inverse(), p2.rho_inverse(), n_simd_elements, + "rho_inverse"); + check_eq(p1.kappa(), p2.kappa(), n_simd_elements, "kappa"); check_eq( - p1.kappa_inverse, + p1.kappa_inverse(), (static_cast< typename specfem::datatype::simd::datatype>( 1.0)) / - p2.kappa, - n_simd_elements); - check_eq(p1.rho_vpinverse, - Kokkos::sqrt(p2.rho_inverse * p2.kappa_inverse), - n_simd_elements); + p2.kappa(), + n_simd_elements, "kappa_inverse"); + check_eq( + p1.rho_vpinverse(), + Kokkos::sqrt(p2.rho_inverse() * + (static_cast::datatype>(1.0) / + p2.kappa())), + n_simd_elements, "rho_vpinverse"); } template ; + constexpr int simd_size = specfem::datatype::simd::size(); @@ -379,14 +394,19 @@ void check_to_value(const specfem::compute::properties properties, get_point_property(ielement + j, iz, ix, properties); const type_real value = values_to_store(i); - if (point_property != value) { - std::ostringstream message; + for (int l = 0; l < PointType::nprops; l++) { + if (point_property.data[l] != value) { + std::ostringstream message; - message << "\n \t Error at ispec = " << ielement + j - << ", iz = " << iz << ", ix = " << ix; - message << get_error_message(point_property, value); + message << "\n \t Error in function check_to_value"; - throw std::runtime_error(message.str()); + message << "\n \t Error at ispec = " << ielement + j + << ", iz = " << iz << ", ix = " << ix + << ", iprop = " << l; + message << get_error_message(point_property, value); + + throw std::runtime_error(message.str()); + } } } } @@ -411,6 +431,9 @@ void check_compute_to_mesh( const int ngllz = properties.ngllz; std::vector elements; + using PointType = specfem::point::properties; + for (int ispec = 0; ispec < nspec; ispec++) { if ((element_types.get_medium_tag(ispec) == MediumTag) && (element_types.get_property_tag(ispec) == PropertyTag)) { @@ -445,17 +468,31 @@ void check_compute_to_mesh( std::get >( materials[ispec_mesh]); auto value = material.get_properties(); - if (point_property != value) { - std::ostringstream message; + for (int l = 0; l < PointType::nprops; l++) { + if (point_property.data[l] != value.data[l]) { + std::ostringstream message; - message << "\n \t Error at ispec = " << ielement << ", iz = " << iz - << ", ix = " << ix; + message << "\n \t Error in function check_compute_to_mesh"; - message << get_error_message(value, 0.0, 1); - message << get_error_message(point_property, 0.0, 2); + message << "\n \t Error at ispec = " << ielement << ", iz = " << iz + << ", ix = " << ix << ", iprop = " << l; + message << get_error_message(value, 0.0, 1); + message << get_error_message(point_property, 0.0, 2); - throw std::runtime_error(message.str()); + throw std::runtime_error(message.str()); + } } + // if (point_property != value) { + // std::ostringstream message; + + // message << "\n \t Error at ispec = " << ielement << ", iz = " << iz + // << ", ix = " << ix; + + // message << get_error_message(value, 0.0, 1); + // message << get_error_message(point_property, 0.0, 2); + + // throw std::runtime_error(message.str()); + // } } } } @@ -515,11 +552,14 @@ void check_store_on_host(specfem::compute::properties &properties, const auto index = get_index(ielement, n_simd_elements, iz, ix); const type_real value = values_to_store_h(i); - PointType point(value); + PointType point; + for (int l = 0; l < PointType::nprops; l++) { + point.data[l] = value; + } PointType point_loaded; specfem::compute::store_on_host(index, properties, point); specfem::compute::load_on_host(index, properties, point_loaded); - check_point_properties(point, point_loaded, + check_point_properties(point_loaded, point, n_simd_elements); } } @@ -614,28 +654,33 @@ void check_load_on_device(specfem::compute::properties &properties, for (int lane = 0; lane < n_simd_elements; lane++) { const auto point_property_l = get_point_property(lane, point_property); - if (point_property_l != value_l) { - std::ostringstream message; + for (int l = 0; l < PointType::nprops; l++) { + if (point_property_l.data[l] != value_l) { + std::ostringstream message; - message << "\n \t Error in function load_on_device"; + message << "\n \t Error in function load_on_device"; - message << "\n \t Error at ispec = " << ielement << ", iz = " << 0 - << ", ix = " << 0; - message << get_error_message(point_property_l, value_l); + message << "\n \t Error at ispec = " << ielement + << ", iz = " << 0 << ", ix = " << 0 + << ", iprop = " << l; + message << get_error_message(point_property_l, value_l); - throw std::runtime_error(message.str()); + throw std::runtime_error(message.str()); + } } } } else if constexpr (!using_simd) { - if (point_property != value_l) { - std::ostringstream message; - message << "\n \t Error in function load_on_device"; + for (int l = 0; l < PointType::nprops; l++) { + if (point_property.data[l] != value_l) { + std::ostringstream message; + message << "\n \t Error in function load_on_device"; - message << "\n \t Error at ispec = " << ielement << ", iz = " << 0 - << ", ix = " << 0; - message << get_error_message(point_property, value_l); + message << "\n \t Error at ispec = " << ielement << ", iz = " << 0 + << ", ix = " << 0 << ", iprop = " << l; + message << get_error_message(point_property, value_l); - throw std::runtime_error(message.str()); + throw std::runtime_error(message.str()); + } } } } From e778c3bb932dd31b1650c1c7abe88a7376ec0661 Mon Sep 17 00:00:00 2001 From: Congyue Cui Date: Tue, 25 Feb 2025 13:32:42 -0500 Subject: [PATCH 2/4] Compute derivated perperties on the fly. --- include/medium/properties_container.hpp | 4 - include/point/properties.hpp | 111 +++++++++--------------- 2 files changed, 41 insertions(+), 74 deletions(-) diff --git a/include/medium/properties_container.hpp b/include/medium/properties_container.hpp index 3fc831ce6..89ec9bb99 100644 --- a/include/medium/properties_container.hpp +++ b/include/medium/properties_container.hpp @@ -61,8 +61,6 @@ struct impl_properties_container { for (int i = 0; i < nprops; i++) { property.data[i] = target(ispec, iz, ix, i); } - - property.compute(); } template KOKKOS_FORCEINLINE_FUNCTION void @@ -90,8 +88,6 @@ struct impl_properties_container { Kokkos::Experimental::where(mask, property.data[i]) .copy_from(&target(ispec, iz, ix, i), tag_type()); } - - property.compute(); } public: diff --git a/include/point/properties.hpp b/include/point/properties.hpp index fab5902dd..6b703ddc7 100644 --- a/include/point/properties.hpp +++ b/include/point/properties.hpp @@ -8,22 +8,21 @@ namespace point { #define DEFINE_PROP(prop) \ constexpr static int i_##prop = __COUNTER__ - _counter - 1; \ - KOKKOS_INLINE_FUNCTION value_type prop() const { \ + KOKKOS_INLINE_FUNCTION constexpr value_type prop() const { \ return Base::data[i_##prop]; \ } \ - KOKKOS_INLINE_FUNCTION void prop(value_type val) { \ + KOKKOS_INLINE_FUNCTION constexpr void prop(value_type &val) { \ Base::data[i_##prop] = val; \ } namespace impl { -template struct impl_properties { +template struct impl_properties { using simd = specfem::datatype::simd; ///< SIMD type using value_type = typename simd::datatype; constexpr static bool is_point_properties = true; constexpr static auto nprops = N; - constexpr static auto nprops_extra = N_EX; - value_type data[N_EX]; + value_type data[N]; KOKKOS_FUNCTION impl_properties() = default; @@ -52,8 +51,8 @@ template struct impl_properties { * */ KOKKOS_FUNCTION - bool operator==(const impl_properties &rhs) const { - for (int i = 0; i < N_EX; ++i) { + bool operator==(const impl_properties &rhs) const { + for (int i = 0; i < N; ++i) { if (data[i] != rhs.data[i]) { return false; } @@ -66,7 +65,7 @@ template struct impl_properties { * */ KOKKOS_FUNCTION - bool operator!=(const impl_properties &rhs) const { + bool operator!=(const impl_properties &rhs) const { return !(*this == rhs); } }; @@ -94,14 +93,13 @@ template struct properties - : public impl::impl_properties<3, 6, UseSIMD> { + : public impl::impl_properties<3, UseSIMD> { /** * @name Typedefs * */ ///@{ - using Base = impl::impl_properties<3, 6, UseSIMD>; constexpr static auto dimension = specfem::dimension::type::dim2; constexpr static auto medium_tag = specfem::element::medium_tag::elastic; constexpr static auto property_tag = @@ -111,31 +109,25 @@ struct properties = 0> - KOKKOS_FUNCTION properties(Args... args) : Base(args...) { - compute(); - } + using Base = impl::impl_properties<3, UseSIMD>; + using Base::Base; DEFINE_PROP(lambdaplus2mu) ///< Lame's parameter @f$ \lambda + 2\mu @f$ DEFINE_PROP(mu) ///< shear modulus @f$ \mu @f$ DEFINE_PROP(rho) ///< density @f$ \rho @f$ - DEFINE_PROP(rho_vp) ///< P-wave velocity @f$ \rho v_p @f$ - DEFINE_PROP(rho_vs) ///< S-wave velocity @f$ \rho v_s @f$ - DEFINE_PROP(lambda) ///< Lame's parameter @f$ \lambda @f$ + KOKKOS_INLINE_FUNCTION constexpr value_type rho_vp() const { + return Kokkos::sqrt(rho() * lambdaplus2mu()); ///< P-wave velocity @f$ \rho + ///< v_p @f$ + } - KOKKOS_INLINE_FUNCTION - void compute() { - rho_vp(Kokkos::sqrt(rho() * lambdaplus2mu())); - rho_vs(Kokkos::sqrt(rho() * mu())); - lambda(lambdaplus2mu() - (static_cast(2.0)) * mu()); + KOKKOS_INLINE_FUNCTION constexpr value_type rho_vs() const { + return Kokkos::sqrt(rho() * mu()); ///< S-wave velocity @f$ \rho v_s @f$ + } + + KOKKOS_INLINE_FUNCTION constexpr value_type lambda() const { + return lambdaplus2mu() - (static_cast(2.0)) * + mu(); ///< Lame's parameter @f$ \lambda @f$ } }; @@ -143,14 +135,13 @@ template struct properties - : public impl::impl_properties<10, 12, UseSIMD> { + : public impl::impl_properties<10, UseSIMD> { /** * @name Typedefs * */ ///@{ - using Base = impl::impl_properties<10, 12, UseSIMD>; constexpr static auto dimension = specfem::dimension::type::dim2; constexpr static auto medium_tag = specfem::element::medium_tag::elastic; constexpr static auto property_tag = @@ -160,17 +151,8 @@ struct properties = 0> - KOKKOS_FUNCTION properties(Args... args) : Base(args...) { - compute(); - } + using Base = impl::impl_properties<10, UseSIMD>; + using Base::Base; /** * @name Elastic constants @@ -186,16 +168,15 @@ struct properties struct properties - : public impl::impl_properties<2, 4, UseSIMD> { - + : public impl::impl_properties<2, UseSIMD> { /** * @name Typedefs * */ ///@{ - using Base = impl::impl_properties<2, 4, UseSIMD>; constexpr static auto dimension = specfem::dimension::type::dim2; constexpr static auto medium_tag = specfem::element::medium_tag::acoustic; constexpr static auto property_tag = @@ -225,31 +204,23 @@ struct properties; + using Base::Base; - KOKKOS_FUNCTION - properties(const value_type *value) : Base(value) { compute(); } + DEFINE_PROP(rho_inverse) ///< @f$ \frac{1}{\rho} @f$ + DEFINE_PROP(kappa) ///< Bulk modulus @f$ \kappa @f$ - template = 0> - KOKKOS_FUNCTION properties(Args... args) : Base(args...) { - compute(); + KOKKOS_INLINE_FUNCTION constexpr value_type kappa_inverse() const { + return (static_cast(1.0)) / + kappa(); ///< @f$ \frac{1}{\lambda + 2\mu} @f$ } - DEFINE_PROP(rho_inverse) ///< @f$ \frac{1}{\rho} @f$ - DEFINE_PROP(kappa) ///< Bulk modulus @f$ \kappa @f$ - DEFINE_PROP(kappa_inverse) ///< @f$ \frac{1}{\lambda + 2\mu} @f$ - DEFINE_PROP(rho_vpinverse) ///< @f$ \frac{1}{\rho v_p} @f$ - - KOKKOS_INLINE_FUNCTION - void compute() { - kappa_inverse((static_cast(1.0)) / kappa()); - rho_vpinverse(Kokkos::sqrt(rho_inverse() * kappa_inverse())); + KOKKOS_INLINE_FUNCTION constexpr value_type rho_vpinverse() const { + return Kokkos::sqrt(rho_inverse() * kappa_inverse()); ///< @f$ \frac{1}{\rho + ///< v_p} @f$ } }; -#undef DEFINE_PROP_CONSTRUCTORS #undef DEFINE_PROP } // namespace point From c09796f51858d29bc9d34a20aeb98edc5eda9872 Mon Sep 17 00:00:00 2001 From: Congyue Cui Date: Tue, 25 Feb 2025 14:46:50 -0500 Subject: [PATCH 3/4] Update macro location. --- include/point/properties.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/point/properties.hpp b/include/point/properties.hpp index 6b703ddc7..dfa3c9015 100644 --- a/include/point/properties.hpp +++ b/include/point/properties.hpp @@ -3,9 +3,6 @@ #include "datatypes/simd.hpp" #include "enumerations/medium.hpp" -namespace specfem { -namespace point { - #define DEFINE_PROP(prop) \ constexpr static int i_##prop = __COUNTER__ - _counter - 1; \ KOKKOS_INLINE_FUNCTION constexpr value_type prop() const { \ @@ -15,6 +12,9 @@ namespace point { Base::data[i_##prop] = val; \ } +namespace specfem { +namespace point { + namespace impl { template struct impl_properties { using simd = specfem::datatype::simd; ///< SIMD type @@ -221,7 +221,7 @@ struct properties Date: Wed, 26 Feb 2025 10:18:16 -0500 Subject: [PATCH 4/4] Rename Base to base_type --- .../isotropic/properties_container.hpp | 4 +-- .../anisotropic/properties_container.hpp | 4 +-- .../isotropic/properties_container.hpp | 4 +-- include/medium/properties_container.hpp | 4 +-- include/point/properties.hpp | 16 ++++----- include/policies/chunk.hpp | 34 +++++++++---------- 6 files changed, 33 insertions(+), 33 deletions(-) diff --git a/include/medium/dim2/acoustic/isotropic/properties_container.hpp b/include/medium/dim2/acoustic/isotropic/properties_container.hpp index ee9381fec..6a7881386 100644 --- a/include/medium/dim2/acoustic/isotropic/properties_container.hpp +++ b/include/medium/dim2/acoustic/isotropic/properties_container.hpp @@ -13,11 +13,11 @@ struct properties_container { - using Base = + using base_type = impl::impl_properties_container; - using Base::Base; + using base_type::base_type; constexpr static int _counter = __COUNTER__; DEFINE_CONTAINER(rho_inverse) diff --git a/include/medium/dim2/elastic/anisotropic/properties_container.hpp b/include/medium/dim2/elastic/anisotropic/properties_container.hpp index b13079853..b78debd9d 100644 --- a/include/medium/dim2/elastic/anisotropic/properties_container.hpp +++ b/include/medium/dim2/elastic/anisotropic/properties_container.hpp @@ -13,10 +13,10 @@ struct properties_container { - using Base = impl::impl_properties_container< + using base_type = impl::impl_properties_container< specfem::element::medium_tag::elastic, specfem::element::property_tag::anisotropic, 10>; - using Base::Base; + using base_type::base_type; constexpr static int _counter = __COUNTER__; DEFINE_CONTAINER(c11) diff --git a/include/medium/dim2/elastic/isotropic/properties_container.hpp b/include/medium/dim2/elastic/isotropic/properties_container.hpp index 1a4e8bdcf..998d88776 100644 --- a/include/medium/dim2/elastic/isotropic/properties_container.hpp +++ b/include/medium/dim2/elastic/isotropic/properties_container.hpp @@ -13,11 +13,11 @@ struct properties_container { - using Base = + using base_type = impl::impl_properties_container; - using Base::Base; + using base_type::base_type; constexpr static int _counter = __COUNTER__; DEFINE_CONTAINER(lambdaplus2mu) diff --git a/include/medium/properties_container.hpp b/include/medium/properties_container.hpp index 89ec9bb99..9a762252b 100644 --- a/include/medium/properties_container.hpp +++ b/include/medium/properties_container.hpp @@ -6,11 +6,11 @@ constexpr static int i_##prop = __COUNTER__ - _counter - 1; \ KOKKOS_INLINE_FUNCTION type_real &prop(const int &ispec, const int &iz, \ const int &ix) const { \ - return Base::data(ispec, iz, ix, i_##prop); \ + return base_type::data(ispec, iz, ix, i_##prop); \ } \ KOKKOS_INLINE_FUNCTION type_real &h_##prop(const int &ispec, const int &iz, \ const int &ix) const { \ - return Base::h_data(ispec, iz, ix, i_##prop); \ + return base_type::h_data(ispec, iz, ix, i_##prop); \ } namespace specfem { diff --git a/include/point/properties.hpp b/include/point/properties.hpp index dfa3c9015..34befd8b4 100644 --- a/include/point/properties.hpp +++ b/include/point/properties.hpp @@ -6,10 +6,10 @@ #define DEFINE_PROP(prop) \ constexpr static int i_##prop = __COUNTER__ - _counter - 1; \ KOKKOS_INLINE_FUNCTION constexpr value_type prop() const { \ - return Base::data[i_##prop]; \ + return base_type::data[i_##prop]; \ } \ KOKKOS_INLINE_FUNCTION constexpr void prop(value_type &val) { \ - Base::data[i_##prop] = val; \ + base_type::data[i_##prop] = val; \ } namespace specfem { @@ -109,8 +109,8 @@ struct properties; - using Base::Base; + using base_type = impl::impl_properties<3, UseSIMD>; + using base_type::base_type; DEFINE_PROP(lambdaplus2mu) ///< Lame's parameter @f$ \lambda + 2\mu @f$ DEFINE_PROP(mu) ///< shear modulus @f$ \mu @f$ @@ -151,8 +151,8 @@ struct properties; - using Base::Base; + using base_type = impl::impl_properties<10, UseSIMD>; + using base_type::base_type; /** * @name Elastic constants @@ -204,8 +204,8 @@ struct properties; - using Base::Base; + using base_type = impl::impl_properties<2, UseSIMD>; + using base_type::base_type; DEFINE_PROP(rho_inverse) ///< @f$ \frac{1}{\rho} @f$ DEFINE_PROP(kappa) ///< Bulk modulus @f$ \kappa @f$ diff --git a/include/policies/chunk.hpp b/include/policies/chunk.hpp index 018b99a62..a112f5f51 100644 --- a/include/policies/chunk.hpp +++ b/include/policies/chunk.hpp @@ -63,14 +63,14 @@ struct chunk_index_type { template struct mapped_chunk_index_type : public chunk_index_type { - using Base = chunk_index_type; + using base_type = chunk_index_type; int imap; ///< Index of the mapped element KOKKOS_INLINE_FUNCTION mapped_chunk_index_type(const int ielement, const specfem::point::index index, const int imap) - : Base(ielement, index), imap(imap) {} + : base_type(ielement, index), imap(imap) {} }; } // namespace impl @@ -252,20 +252,20 @@ class mapped_chunk; template class mapped_chunk : public chunk { - using Base = chunk; + using base_type = chunk; using mapped_index_type = - typename impl::mapped_chunk_index_type; ///< Index + typename impl::mapped_chunk_index_type; ///< Index public: KOKKOS_INLINE_FUNCTION mapped_chunk(const ViewType &indices, const ViewType &mapping, const int ngllz, const int ngllx) - : Base(indices, ngllz, ngllx), mapping(mapping) {} + : base_type(indices, ngllz, ngllx), mapping(mapping) {} KOKKOS_INLINE_FUNCTION mapped_index_type operator()(const int i) const { - const auto base_index = Base::operator()(i); + const auto base_index = base_type::operator()(i); return mapped_index_type(base_index.ielement, base_index.index, mapping(base_index.ielement)); } @@ -405,8 +405,8 @@ struct element_chunk template struct mapped_element_chunk : public element_chunk { using simd = typename ParallelConfig::simd; - using Base = element_chunk; - using IndexViewType = typename Base::IndexViewType; + using base_type = element_chunk; + using IndexViewType = typename base_type::IndexViewType; using mapped_iterator_type = specfem::iterator::mapped_chunk { mapped_element_chunk(const IndexViewType &view, const IndexViewType &mapping, int ngllz, int ngllx) - : Base(view, ngllz, ngllx), mapping(mapping) {} + : base_type(view, ngllz, ngllx), mapping(mapping) {} KOKKOS_INLINE_FUNCTION mapped_iterator_type mapped_league_iterator(const int start_index) const { const int start = start_index; - const int end = - (start + Base::chunk_size * Base::simd_size > Base::elements.extent(0)) - ? Base::elements.extent(0) - : start + Base::chunk_size * Base::simd_size; + const int end = (start + base_type::chunk_size * base_type::simd_size > + base_type::elements.extent(0)) + ? base_type::elements.extent(0) + : start + base_type::chunk_size * base_type::simd_size; const auto elem_indices = - Kokkos::subview(Base::elements, Kokkos::make_pair(start, end)); + Kokkos::subview(base_type::elements, Kokkos::make_pair(start, end)); const auto map_indices = Kokkos::subview(mapping, Kokkos::make_pair(start, end)); - return mapped_iterator_type(elem_indices, map_indices, Base::ngllz, - Base::ngllx); + return mapped_iterator_type(elem_indices, map_indices, base_type::ngllz, + base_type::ngllx); } protected: