Skip to content

Commit

Permalink
Update point_property implementation
Browse files Browse the repository at this point in the history
  • Loading branch information
@icui
icui committed Feb 19, 2025
1 parent 337b147 commit d83f142
Show file tree
Hide file tree
Showing 3 changed files with 112 additions and 146 deletions.
114 changes: 41 additions & 73 deletions benchmarks/impl_point.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
namespace specfem {
namespace benchmarks {

#define PROPERTIES_CONSTRUCTOR \
#define DEFINE_PROP_CONSTRUCTORS \
KOKKOS_FUNCTION \
properties() = default; \
KOKKOS_FUNCTION \
Expand All @@ -18,6 +18,10 @@ namespace benchmarks {
compute(); \
}

#define DEFINE_PROP(prop, i) \
KOKKOS_INLINE_FUNCTION value_type prop() const { return Base::data[i]; } \
KOKKOS_INLINE_FUNCTION void prop(value_type val) { Base::data[i] = val; }

namespace impl {
template <int N, int NALL, bool UseSIMD> struct impl_properties {
using simd = specfem::datatype::simd<type_real, UseSIMD>; ///< SIMD type
Expand Down Expand Up @@ -114,31 +118,21 @@ struct properties<specfem::dimension::type::dim2,
typename simd::datatype; ///< Value type to store properties
///@}

value_type &lambdaplus2mu =
Base::data[0]; ///< Lame's parameter @f$ \lambda + 2\mu @f$
value_type &mu = Base::data[1]; ///< shear modulus @f$ \mu @f$
value_type &rho = Base::data[2]; ///< density @f$ \rho @f$

value_type &rho_vp = Base::data[3]; ///< P-wave velocity @f$ \rho v_p @f$
value_type &rho_vs = Base::data[4]; ///< S-wave velocity @f$ \rho v_s @f$
value_type &lambda = Base::data[5]; ///< Lame's parameter @f$ \lambda @f$
DEFINE_PROP(lambdaplus2mu, 0) ///< Lame's parameter @f$ \lambda + 2\mu @f$
DEFINE_PROP(mu, 1) ///< shear modulus @f$ \mu @f$
DEFINE_PROP(rho, 2) ///< density @f$ \rho @f$

PROPERTIES_CONSTRUCTOR
DEFINE_PROP(rho_vp, 3) ///< P-wave velocity @f$ \rho v_p @f$
DEFINE_PROP(rho_vs, 4) ///< S-wave velocity @f$ \rho v_s @f$
DEFINE_PROP(lambda, 5) ///< Lame's parameter @f$ \lambda @f$

// KOKKOS_FUNCTION
// properties(const value_type &lambdaplus2mu, const value_type &mu,
// const value_type &rho) {
// this->lambdaplus2mu = lambdaplus2mu;
// this->mu = mu;
// this->rho = rho;
// compute();
// }
DEFINE_PROP_CONSTRUCTORS

KOKKOS_INLINE_FUNCTION
void compute() {
rho_vp = Kokkos::sqrt(rho * lambdaplus2mu);
rho_vs = Kokkos::sqrt(rho * mu);
lambda = lambdaplus2mu - (static_cast<value_type>(2.0)) * mu;
rho_vp(Kokkos::sqrt(rho() * lambdaplus2mu()));
rho_vs(Kokkos::sqrt(rho() * mu()));
lambda(lambdaplus2mu() - (static_cast<value_type>(2.0)) * mu());
}
};

Expand Down Expand Up @@ -169,46 +163,27 @@ struct properties<specfem::dimension::type::dim2,
*
*/
///@{
value_type &c11 = Base::data[0]; ///< @f$ c_{11} @f$
value_type &c13 = Base::data[1]; ///< @f$ c_{13} @f$
value_type &c15 = Base::data[2]; ///< @f$ c_{15} @f$
value_type &c33 = Base::data[3]; ///< @f$ c_{33} @f$
value_type &c35 = Base::data[4]; ///< @f$ c_{35} @f$
value_type &c55 = Base::data[5]; ///< @f$ c_{55} @f$
value_type &c12 = Base::data[6]; ///< @f$ c_{12} @f$
value_type &c23 = Base::data[7]; ///< @f$ c_{23} @f$
value_type &c25 = Base::data[8]; ///< @f$ c_{25} @f$
DEFINE_PROP(c11, 0); ///< @f$ c_{11} @f$
DEFINE_PROP(c13, 1); ///< @f$ c_{13} @f$
DEFINE_PROP(c15, 2); ///< @f$ c_{15} @f$
DEFINE_PROP(c33, 3); ///< @f$ c_{33} @f$
DEFINE_PROP(c35, 4); ///< @f$ c_{35} @f$
DEFINE_PROP(c55, 5); ///< @f$ c_{55} @f$
DEFINE_PROP(c12, 6); ///< @f$ c_{12} @f$
DEFINE_PROP(c23, 7); ///< @f$ c_{23} @f$
DEFINE_PROP(c25, 8); ///< @f$ c_{25} @f$
///@}

value_type &rho = Base::data[9]; ///< Density @f$ \rho @f$
value_type &rho_vp = Base::data[10]; ///< P-wave velocity @f$ \rho v_p @f$
value_type &rho_vs = Base::data[11]; ///< S-wave velocity @f$ \rho v_s @f$

PROPERTIES_CONSTRUCTOR

// 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) {
// this->c11 = c11;
// this->c13 = c13;
// this->c15 = c15;
// this->c33 = c33;
// this->c35 = c35;
// this->c55 = c55;
// this->c12 = c12;
// this->c23 = c23;
// this->c25 = c25;
// this->rho = rho;
// compute();
// }
DEFINE_PROP(rho, 9); ///< Density @f$ \rho @f$
DEFINE_PROP(rho_vp, 10); ///< P-wave velocity @f$ \rho v_p @f$
DEFINE_PROP(rho_vs, 11); ///< S-wave velocity @f$ \rho v_s @f$

DEFINE_PROP_CONSTRUCTORS

KOKKOS_INLINE_FUNCTION
void compute() {
rho_vp = Kokkos::sqrt(rho * c33);
rho_vs = Kokkos::sqrt(rho * c55);
rho_vp(Kokkos::sqrt(rho() * c33()));
rho_vs(Kokkos::sqrt(rho() * c55()));
}
};

Expand Down Expand Up @@ -240,30 +215,23 @@ struct properties<specfem::dimension::type::dim2,
typename simd::datatype; ///< Value type to store properties
///@}

value_type &rho_inverse = Base::data[0]; ///< @f$ \frac{1}{\rho} @f$
value_type &kappa = Base::data[1]; ///< Bulk modulus @f$ \kappa @f$

value_type &kappa_inverse =
Base::data[2]; ///< @f$ \frac{1}{\lambda + 2\mu} @f$
value_type &rho_vpinverse = Base::data[3]; ///< @f$ \frac{1}{\rho v_p} @f$

PROPERTIES_CONSTRUCTOR
DEFINE_PROP_CONSTRUCTORS

// KOKKOS_FUNCTION
// properties(const value_type &rho_inverse, const value_type &kappa) {
// this->rho_inverse = rho_inverse;
// this->kappa = kappa;
// compute();
// }
DEFINE_PROP(rho_inverse, 0) ///< @f$ \frac{1}{\rho} @f$
DEFINE_PROP(kappa, 1) ///< Bulk modulus @f$ \kappa @f$
DEFINE_PROP(kappa_inverse, 2) ///< @f$ \frac{1}{\lambda + 2\mu} @f$
DEFINE_PROP(rho_vpinverse, 3) ///< @f$ \frac{1}{\rho v_p} @f$

private:
KOKKOS_INLINE_FUNCTION
void compute() {
kappa_inverse = (static_cast<value_type>(1.0)) / kappa;
rho_vpinverse = Kokkos::sqrt(rho_inverse * kappa_inverse);
kappa_inverse((static_cast<value_type>(1.0)) / kappa());
rho_vpinverse(Kokkos::sqrt(rho_inverse() * kappa_inverse()));
}
};

#undef PROPERTIES_CONSTRUCTOR
#undef DEFINE_PROP_CONSTRUCTORS
#undef DEFINE_PROP

} // namespace benchmarks
} // namespace specfem
34 changes: 19 additions & 15 deletions benchmarks/main.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// #include "archive/stiffness1.hpp"
#include "archive/stiffness1.hpp"
#include "execute.hpp"
#include "stiffness.hpp"
// #include "divide.hpp"
Expand Down Expand Up @@ -47,26 +47,30 @@ void benchmark(specfem::compute::assembly &assembly,
// field.elastic.nglob, field.elastic.components);

for (const auto [istep, dt] : time_scheme->iterate_forward()) {
compute_stiffness_interaction<acoustic, isotropic, flag>(assembly, istep);
compute_stiffness_interaction<elastic, isotropic, flag>(assembly, istep);
compute_stiffness_interaction<elastic, anisotropic, flag>(assembly, istep);
// compute_stiffness_interaction<acoustic, isotropic, flag>(assembly,
// istep); compute_stiffness_interaction<elastic, isotropic, flag>(assembly,
// istep); compute_stiffness_interaction<elastic, anisotropic,
// flag>(assembly, istep);

// compute_stiffness_interaction2<acoustic, isotropic, flag>(assembly,
// istep); compute_stiffness_interaction2<elastic, isotropic,
// flag>(assembly, istep); compute_stiffness_interaction2<elastic,
// anisotropic, flag>(assembly, istep);

// if constexpr (flag) {
// compute_stiffness_interaction<acoustic, isotropic, false>(assembly,
// istep); compute_stiffness_interaction<elastic, isotropic,
// false>(assembly, istep); compute_stiffness_interaction<elastic,
// anisotropic, false>(assembly, istep);
// } else {
// compute_stiffness_interaction2<acoustic, isotropic, false>(assembly,
// istep); compute_stiffness_interaction2<elastic, isotropic,
// false>(assembly, istep); compute_stiffness_interaction2<elastic,
// anisotropic, false>(assembly, istep);
// }
if constexpr (flag) {
compute_stiffness_interaction<acoustic, isotropic, false>(assembly,
istep);
compute_stiffness_interaction<elastic, isotropic, false>(assembly, istep);
compute_stiffness_interaction<elastic, anisotropic, false>(assembly,
istep);
} else {
compute_stiffness_interaction2<acoustic, isotropic, false>(assembly,
istep);
compute_stiffness_interaction2<elastic, isotropic, false>(assembly,
istep);
compute_stiffness_interaction2<elastic, anisotropic, false>(assembly,
istep);
}

// divide_mass_matrix<dimension, wavefield, acoustic>(assembly);
// divide_mass_matrix<dimension, wavefield, elastic>(assembly);
Expand Down
110 changes: 52 additions & 58 deletions benchmarks/stiffness.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,15 +187,9 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
xiz * df_dxi[icomponent] + gammaz * df_dgamma[icomponent];
}

auto point_property = [&]() {
if constexpr (flag) {
return specfem::point::properties<
dimension, medium_tag, property_tag, using_simd>();
} else {
return specfem::benchmarks::properties<
dimension, medium_tag, property_tag, using_simd>();
}
}();
specfem::benchmarks::properties<dimension, medium_tag,
property_tag, using_simd>
point_property;

specfem::point::field_derivatives<dimension, medium_tag,
using_simd>
Expand All @@ -211,23 +205,23 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
specfem::element::medium_tag::acoustic &&
PropertyTag ==
specfem::element::property_tag::isotropic) {
point_property.rho_inverse =
container.rho_inverse(ispec, iz, ix);
point_property.kappa = container.kappa(ispec, iz, ix);
point_property.kappa_inverse =
static_cast<type_real>(1.0) / point_property.kappa;
point_property.rho_vpinverse =
sqrt(point_property.rho_inverse *
point_property.kappa_inverse);
point_property.rho_inverse(
container.rho_inverse(ispec, iz, ix));
point_property.kappa(container.kappa(ispec, iz, ix));
point_property.kappa_inverse(static_cast<type_real>(1.0) /
point_property.kappa());
point_property.rho_vpinverse(
sqrt(point_property.rho_inverse() *
point_property.kappa_inverse()));

const auto &du = field_derivatives.du;

specfem::datatype::VectorPointViewType<type_real, 2, 1,
using_simd>
T;

T(0, 0) = point_property.rho_inverse * du(0, 0);
T(1, 0) = point_property.rho_inverse * du(1, 0);
T(0, 0) = point_property.rho_inverse() * du(0, 0);
T(1, 0) = point_property.rho_inverse() * du(1, 0);

return specfem::point::stress<
specfem::dimension::type::dim2,
Expand All @@ -237,16 +231,16 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
MediumTag == specfem::element::medium_tag::elastic &&
PropertyTag ==
specfem::element::property_tag::isotropic) {
point_property.rho = container.rho(ispec, iz, ix);
point_property.mu = container.mu(ispec, iz, ix);
point_property.lambdaplus2mu =
container.lambdaplus2mu(ispec, iz, ix);
point_property.lambda =
point_property.lambdaplus2mu - 2 * point_property.mu;
point_property.rho_vp =
sqrt(point_property.rho * point_property.lambdaplus2mu);
point_property.rho_vs =
sqrt(point_property.rho * point_property.mu);
point_property.rho(container.rho(ispec, iz, ix));
point_property.mu(container.mu(ispec, iz, ix));
point_property.lambdaplus2mu(
container.lambdaplus2mu(ispec, iz, ix));
point_property.lambda(point_property.lambdaplus2mu() -
2 * point_property.mu());
point_property.rho_vp(sqrt(point_property.rho() *
point_property.lambdaplus2mu()));
point_property.rho_vs(
sqrt(point_property.rho() * point_property.mu()));

using datatype =
typename specfem::datatype::simd<type_real,
Expand All @@ -255,11 +249,11 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,

datatype sigma_xx, sigma_zz, sigma_xz;

sigma_xx = point_property.lambdaplus2mu * du(0, 0) +
point_property.lambda * du(1, 1);
sigma_zz = point_property.lambdaplus2mu * du(1, 1) +
point_property.lambda * du(0, 0);
sigma_xz = point_property.mu * (du(0, 1) + du(1, 0));
sigma_xx = point_property.lambdaplus2mu() * du(0, 0) +
point_property.lambda() * du(1, 1);
sigma_zz = point_property.lambdaplus2mu() * du(1, 1) +
point_property.lambda() * du(0, 0);
sigma_xz = point_property.mu() * (du(0, 1) + du(1, 0));

specfem::datatype::VectorPointViewType<type_real, 2, 2,
using_simd>
Expand All @@ -277,21 +271,21 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,
MediumTag == specfem::element::medium_tag::elastic &&
PropertyTag ==
specfem::element::property_tag::anisotropic) {
point_property.rho = container.rho(ispec, iz, ix);
point_property.c11 = container.c11(ispec, iz, ix);
point_property.c12 = container.c12(ispec, iz, ix);
point_property.c13 = container.c13(ispec, iz, ix);
point_property.c15 = container.c15(ispec, iz, ix);
point_property.c33 = container.c33(ispec, iz, ix);
point_property.c35 = container.c35(ispec, iz, ix);
point_property.c55 = container.c55(ispec, iz, ix);
point_property.c23 = container.c23(ispec, iz, ix);
point_property.c25 = container.c25(ispec, iz, ix);

point_property.rho_vp =
sqrt(point_property.rho * point_property.c33);
point_property.rho_vs =
sqrt(point_property.rho * point_property.c55);
point_property.rho(container.rho(ispec, iz, ix));
point_property.c11(container.c11(ispec, iz, ix));
point_property.c12(container.c12(ispec, iz, ix));
point_property.c13(container.c13(ispec, iz, ix));
point_property.c15(container.c15(ispec, iz, ix));
point_property.c33(container.c33(ispec, iz, ix));
point_property.c35(container.c35(ispec, iz, ix));
point_property.c55(container.c55(ispec, iz, ix));
point_property.c23(container.c23(ispec, iz, ix));
point_property.c25(container.c25(ispec, iz, ix));

point_property.rho_vp(
sqrt(point_property.rho() * point_property.c33()));
point_property.rho_vs(
sqrt(point_property.rho() * point_property.c55()));

using datatype =
typename specfem::datatype::simd<type_real,
Expand All @@ -300,17 +294,17 @@ void compute_stiffness_interaction(const specfem::compute::assembly &assembly,

datatype sigma_xx, sigma_zz, sigma_xz;

sigma_xx = point_property.c11 * du(0, 0) +
point_property.c13 * du(1, 1) +
point_property.c15 * (du(1, 0) + du(0, 1));
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 = point_property.c13 * du(0, 0) +
point_property.c33 * du(1, 1) +
point_property.c35 * (du(1, 0) + du(0, 1));
sigma_zz = point_property.c13() * du(0, 0) +
point_property.c33() * du(1, 1) +
point_property.c35() * (du(1, 0) + du(0, 1));

sigma_xz = point_property.c15 * du(0, 0) +
point_property.c35 * du(1, 1) +
point_property.c55 * (du(1, 0) + du(0, 1));
sigma_xz = point_property.c15() * du(0, 0) +
point_property.c35() * du(1, 1) +
point_property.c55() * (du(1, 0) + du(0, 1));

specfem::datatype::VectorPointViewType<type_real, 2, 2,
using_simd>
Expand Down

0 comments on commit d83f142

Please sign in to comment.