Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

refactor point_properties and property_container #509

Open
wants to merge 4 commits into
base: issue-493
Choose a base branch
from
Open

refactor point_properties and property_container #509

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
240222e
refactor point_properties and property_container
icui Feb 24, 2025
e778c3b
Compute derivated perperties on the fly.
icui Feb 25, 2025
c09796f
Update macro location.
icui Feb 25, 2025
558545b
Rename Base to base_type
icui Feb 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 3 additions & 15 deletions include/IO/property/reader.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,31 +23,19 @@ void specfem::IO::property_reader<InputLibrary>::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"
Expand Down
18 changes: 3 additions & 15 deletions include/IO/property/writer.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,9 +57,7 @@ void specfem::IO::property_writer<OutputLibrary>::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();
}

{
Expand All @@ -86,16 +84,7 @@ void specfem::IO::property_writer<OutputLibrary>::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();
}

{
Expand All @@ -120,8 +109,7 @@ void specfem::IO::property_writer<OutputLibrary>::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);
Expand Down
20 changes: 10 additions & 10 deletions include/boundary_conditions/stacey/stacey.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ impl_enforce_traction(const acoustic_type &, const isotropic_type &,

// Apply Stacey boundary condition
traction(0) += static_cast<type_real>(-1.0) * factor *
property.rho_vpinverse * field.velocity(0);
property.rho_vpinverse() * field.velocity(0);

// Do nothing
return;
Expand Down Expand Up @@ -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<type_real>(-1.0) * factor *
property.rho_vpinverse * field.velocity(0);
property.rho_vpinverse() * field.velocity(0);

return;
}
Expand Down Expand Up @@ -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<type_real>(-1.0) * factor[0] * jacobian1d *
Expand Down Expand Up @@ -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)) =
Expand Down Expand Up @@ -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<type_real>(-1.0) * factor[0] * jacobian1d *
Expand Down Expand Up @@ -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)) =
Expand Down
4 changes: 2 additions & 2 deletions include/medium/dim2/acoustic/isotropic/frechet_derivative.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<type_real>(1.0) / properties.kappa * dt;
static_cast<type_real>(1.0) / properties.kappa() * dt;

return { rho_kl, kappa_kl };
}
Expand Down
2 changes: 1 addition & 1 deletion include/medium/dim2/acoustic/isotropic/mass_matrix.tpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,5 +14,5 @@ specfem::medium::impl_mass_matrix_component(
specfem::dimension::type::dim2, true, UseSIMD> &partial_derivatives) {

return specfem::datatype::ScalarPointViewType<type_real, 1, UseSIMD>(
partial_derivatives.jacobian / properties.kappa);
partial_derivatives.jacobian / properties.kappa());
}
218 changes: 13 additions & 205 deletions include/medium/dim2/acoustic/isotropic/properties_container.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,211 +9,19 @@ namespace medium {

template <>
struct properties_container<specfem::element::medium_tag::acoustic,
specfem::element::property_tag::isotropic> {

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<type_real ***, Kokkos::LayoutLeft,
Kokkos::DefaultExecutionSpace>;

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<!PointProperties::simd::using_simd, int> = 0>
KOKKOS_FORCEINLINE_FUNCTION void
load_device_properties(const specfem::point::index<dimension> &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<type_real>(1.0) / property.kappa;
property.rho_vpinverse =
sqrt(property.rho_inverse * property.kappa_inverse);
}

template <
typename PointProperties,
typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
KOKKOS_FORCEINLINE_FUNCTION void
load_device_properties(const specfem::point::simd_index<dimension> &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<type_real>(1.0) / property.kappa;
property.rho_vpinverse =
Kokkos::sqrt(property.rho_inverse * property.kappa_inverse);
}

template <
typename PointProperties,
typename std::enable_if_t<!PointProperties::simd::using_simd, int> = 0>
inline void
load_host_properties(const specfem::point::index<dimension> &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<type_real>(1.0) / property.kappa;
property.rho_vpinverse =
sqrt(property.rho_inverse * property.kappa_inverse);
}

template <
typename PointProperties,
typename std::enable_if_t<PointProperties::simd::using_simd, int> = 0>
inline void
load_host_properties(const specfem::point::simd_index<dimension> &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<type_real>(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<!PointProperties::simd::using_simd, int> = 0>
inline void assign(const specfem::point::index<dimension> &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<PointProperties::simd::using_simd, int> = 0>
inline void assign(const specfem::point::simd_index<dimension> &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 =
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can we use base_type here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

updated

impl::impl_properties_container<specfem::element::medium_tag::acoustic,
specfem::element::property_tag::isotropic,
2>;
using Base::Base;
constexpr static int _counter = __COUNTER__;

DEFINE_CONTAINER(rho_inverse)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is awesome!

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, this is nuts.

DEFINE_CONTAINER(kappa)
};

} // namespace medium
Expand Down
2 changes: 1 addition & 1 deletion include/medium/dim2/acoustic/isotropic/source.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down
Loading