Skip to content

Commit

Permalink
Merged compute_seismogram routine
Browse files Browse the repository at this point in the history
  • Loading branch information
@Rohit-Kakodkar
Rohit-Kakodkar committed Aug 24, 2023
1 parent 645f9a2 commit d30b73a
Show file tree
Hide file tree
Showing 7 changed files with 216 additions and 71 deletions.
185 changes: 166 additions & 19 deletions include/domain/domain.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ void initialize_views(
constexpr int components = medium::components;

Kokkos::parallel_for(
"specfem::Domain::acoustic::initiaze_views",
"specfem::domain::domain::initiaze_views",
specfem::kokkos::DeviceMDrange<2, Kokkos::Iterate::Left>(
{ 0, 0 }, { nglob, components }),
KOKKOS_LAMBDA(const int iglob, const int idim) {
Expand Down Expand Up @@ -157,14 +157,14 @@ void assign_elemental_properties(
specfem::kokkos::HostMirror1d<int> h_ispec_domain =
Kokkos::create_mirror_view(ispec_domain);

elements = specfem::kokkos::DeviceView1d<
acoustic_tmp::element_container<acoustic_tmp::element_type<qp_type> > >(
"specfem::domain::domain::elements", nelem_domain);
elements =
specfem::kokkos::DeviceView1d<element_container<element_type<qp_type> > >(
"specfem::domain::domain::elements", nelem_domain);

// Get ispec for each element in this domain
int index = 0;
for (int ispec = 0; ispec < nspec; ispec++) {
if (properties.h_ispec_type(ispec) == specfem::enums::element::acoustic) {
if (properties.h_ispec_type(ispec) == value) {
h_ispec_domain(index) = ispec;
index++;
}
Expand Down Expand Up @@ -329,8 +329,7 @@ void initialize_receivers(
// Check if the receiver is in the domain
int index = 0;
for (int irec = 0; irec < ispec_array.extent(0); irec++) {
if (h_ispec_type(h_ispec_array(irec)) ==
specfem::enums::element::acoustic) {
if (h_ispec_type(h_ispec_array(irec)) == value) {
h_my_receivers(index) = irec;
index++;
}
Expand Down Expand Up @@ -402,19 +401,17 @@ specfem::domain::domain<medium, qp_type>::domain(
specfem::quadrature::quadrature *quadx,
specfem::quadrature::quadrature *quadz)
: field(specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
"specfem::Domain::acoustic::field", nglob,
specfem::enums::element::medium::acoustic::components)),
"specfem::domain::domain::field", nglob, medium::components)),
field_dot(specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
"specfem::Domain::acoustic::field_dot", nglob,
specfem::enums::element::medium::acoustic::components)),
"specfem::domain::domain::field_dot", nglob, medium::components)),
field_dot_dot(
specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
"specfem::Domain::acoustic::field_dot_dot", nglob,
specfem::enums::element::medium::acoustic::components)),
"specfem::domain::domain::field_dot_dot", nglob,
medium::components)),
rmass_inverse(
specfem::kokkos::DeviceView2d<type_real, Kokkos::LayoutLeft>(
"specfem::Domain::acoustic::rmass_inverse", nglob,
specfem::enums::element::medium::acoustic::components)),
"specfem::domain::domain::rmass_inverse", nglob,
medium::components)),
quadrature_points(quadrature_points), compute(compute), quadx(quadx),
quadz(quadz) {

Expand Down Expand Up @@ -477,7 +474,7 @@ void specfem::domain::domain<medium, qp_type>::sync_field(
}

template <class medium, class qp_type>
void specfem::domain::domain<specfem::enums::element::medium::acoustic,
void specfem::domain::domain<medium,
qp_type>::sync_field_dot(specfem::sync::kind
kind) {

Expand Down Expand Up @@ -529,7 +526,7 @@ void specfem::domain::domain<medium, qp_type>::divide_mass_matrix() {
const int nglob = this->rmass_inverse.extent(0);

Kokkos::parallel_for(
"specfem::Domain::acoustic::divide_mass_matrix",
"specfem::domain::domain::divide_mass_matrix",
specfem::kokkos::DeviceRange(0, components * nglob),
KOKKOS_CLASS_LAMBDA(const int in) {
const int iglob = in % nglob;
Expand Down Expand Up @@ -576,7 +573,7 @@ void specfem::domain::domain<medium, qp_type>::compute_stiffness_interaction() {
specfem::enums::axes::z>();

Kokkos::parallel_for(
"specfem::Domain::acoustic::compute_gradients",
"specfem::domain::domain::compute_gradients",
specfem::kokkos::DeviceTeam(this->nelem_domain, NTHREADS, NLANES)
.set_scratch_size(0, Kokkos::PerTeam(scratch_size)),
KOKKOS_CLASS_LAMBDA(
Expand Down Expand Up @@ -734,7 +731,7 @@ void specfem::domain::domain<medium, qp_type>::compute_source_interaction(
const auto ibool = this->compute->ibool;

Kokkos::parallel_for(
"specfem::Domain::acoustic::compute_source_interaction",
"specfem::domain::domain::compute_source_interaction",
specfem::kokkos::DeviceTeam(nsources, Kokkos::AUTO, 1),
KOKKOS_CLASS_LAMBDA(
const specfem::kokkos::DeviceTeam::member_type &team_member) {
Expand Down Expand Up @@ -778,4 +775,154 @@ void specfem::domain::domain<medium, qp_type>::compute_source_interaction(
return;
}

template <typename medium, typename qp_type>
void specfem::domain::domain<medium, qp_type>::compute_seismogram(
const int isig_step) {

// Allocate scratch views for field, field_dot & field_dot_dot. Incase of
// acostic domains when calculating displacement, velocity and acceleration
// seismograms we need to compute the derivatives of the field variables. This
// requires summing over all lagrange derivatives at all quadrature points
// within the element. Scratch views speed up this computation by limiting
// global memory accesses.

constexpr int components = medium::components;
const auto ibool = this->compute->ibool;
const auto hprime_xx = this->quadx->get_hprime();
const auto hprime_zz = this->quadz->get_hprime();
// hprime_xx
int scratch_size = quadrature_points.template shmem_size<
type_real, 1, specfem::enums::axes::x, specfem::enums::axes::x>();

// hprime_zz
scratch_size += quadrature_points.template shmem_size<
type_real, 1, specfem::enums::axes::z, specfem::enums::axes::z>();

// field, field_dot, field_dot_dot
scratch_size +=
3 *
quadrature_points
.template shmem_size<type_real, components, specfem::enums::axes::z,
specfem::enums::axes::x>();

Kokkos::parallel_for(
"specfem::domain::domain::compute_seismogram",
specfem::kokkos::DeviceTeam(this->receivers.extent(0), Kokkos::AUTO, 1)
.set_scratch_size(0, Kokkos::PerTeam(scratch_size)),
KOKKOS_CLASS_LAMBDA(
const specfem::kokkos::DeviceTeam::member_type &team_member) {
int ngllx, ngllz;
quadrature_points.get_ngll(&ngllx, &ngllz);
const int irec = team_member.league_rank();
const auto receiver = this->receivers(irec);
const int ispec = receiver.get_ispec();
const auto seismogram_type = receiver.get_seismogram_type();

// Instantiate shared views
// ----------------------------------------------------------------
auto s_hprime_xx = quadrature_points.template ScratchView<
type_real, 1, specfem::enums::axes::x, specfem::enums::axes::x>(
team_member.team_scratch(0));

auto s_hprime_zz = quadrature_points.template ScratchView<
type_real, 1, specfem::enums::axes::z, specfem::enums::axes::z>(
team_member.team_scratch(0));

auto s_field =
quadrature_points.template ScratchView<type_real, components,
specfem::enums::axes::z,
specfem::enums::axes::x>(
team_member.team_scratch(0));

auto s_field_dot =
quadrature_points.template ScratchView<type_real, components,
specfem::enums::axes::z,
specfem::enums::axes::x>(
team_member.team_scratch(0));

auto s_field_dot_dot =
quadrature_points.template ScratchView<type_real, components,
specfem::enums::axes::z,
specfem::enums::axes::x>(
team_member.team_scratch(0));

// Allocate shared views
// ----------------------------------------------------------------
Kokkos::parallel_for(
quadrature_points.template TeamThreadRange<specfem::enums::axes::x,
specfem::enums::axes::x>(
team_member),
[=](const int xz) {
int iz, ix;
sub2ind(xz, ngllx, iz, ix);
s_hprime_xx(iz, ix) = hprime_xx(iz, ix);
});

Kokkos::parallel_for(
quadrature_points.template TeamThreadRange<specfem::enums::axes::z,
specfem::enums::axes::z>(
team_member),
[=](const int xz) {
int iz, ix;
sub2ind(xz, ngllz, iz, ix);
s_hprime_zz(iz, ix) = hprime_zz(iz, ix);
});

Kokkos::parallel_for(
quadrature_points.template TeamThreadRange<specfem::enums::axes::z,
specfem::enums::axes::x>(
team_member),
[=](const int xz) {
int iz, ix;
sub2ind(xz, ngllx, iz, ix);
int iglob = ibool(ispec, iz, ix);
#ifdef KOKKOS_ENABLE_CUDA
#pragma unroll
#endif
for (int icomponent = 0; icomponent < components; icomponent++) {
s_field(iz, ix, icomponent) = field(iglob, icomponent);
s_field_dot(iz, ix, icomponent) = field_dot(iglob, icomponent);
s_field_dot_dot(iz, ix, icomponent) =
field_dot_dot(iglob, icomponent);
}
});

// Get seismogram field
// ----------------------------------------------------------------

Kokkos::parallel_for(
quadrature_points.template TeamThreadRange<specfem::enums::axes::z,
specfem::enums::axes::x>(
team_member),
[=](const int xz) {
receiver.get_field(xz, isig_step, s_field, s_field_dot,
s_field_dot_dot, s_hprime_xx, s_hprime_zz);
});

// compute seismograms components
//-------------------------------------------------------------------
switch (seismogram_type) {
case specfem::enums::seismogram::type::displacement:
case specfem::enums::seismogram::type::velocity:
case specfem::enums::seismogram::type::acceleration:
dimension::array_type<type_real> seismogram_components;
Kokkos::parallel_reduce(
quadrature_points.template TeamThreadRange<
specfem::enums::axes::z, specfem::enums::axes::x>(
team_member),
[=](const int xz,
dimension::array_type<type_real> &l_seismogram_components) {
receiver.compute_seismogram_components(xz, isig_step,
l_seismogram_components);
},
specfem::kokkos::Sum<dimension::array_type<type_real> >(
seismogram_components));
Kokkos::single(Kokkos::PerTeam(team_member), [=] {
receiver.compute_seismogram(isig_step, seismogram_components);
});
break;
}
});
}

#endif
14 changes: 7 additions & 7 deletions include/domain/impl/receivers/acoustic/acoustic2d.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,17 +19,17 @@ class receiver<specfem::enums::element::dimension::dim2,
using medium = specfem::enums::element::medium::acoustic;
using quadrature_points = qp_type;

template <typename T>
template <typename T, int N>
using ScratchViewType =
typename quadrature_points::template ScratchViewType<T>;
typename quadrature_points::template ScratchViewType<T, N>;

KOKKOS_INLINE_FUNCTION virtual void
get_field(const int xz, const int isig_step,
const ScratchViewType<type_real> field,
const ScratchViewType<type_real> field_dot,
const ScratchViewType<type_real> field_dot_dot,
const ScratchViewType<type_real> hprime_xx,
const ScratchViewType<type_real> hprime_zz) const {};
const ScratchViewType<type_real, medium::components> field,
const ScratchViewType<type_real, medium::components> field_dot,
const ScratchViewType<type_real, medium::components> field_dot_dot,
const ScratchViewType<type_real, 1> hprime_xx,
const ScratchViewType<type_real, 1> hprime_zz) const {};
KOKKOS_INLINE_FUNCTION virtual void compute_seismogram_components(
const int xz, const int isig_step,
dimension::array_type<type_real> &l_seismogram_components) const {};
Expand Down
13 changes: 7 additions & 6 deletions include/domain/impl/receivers/acoustic/acoustic2d_isotropic.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,12 +64,13 @@ class receiver<
sv_receiver_field_type receiver_field);

KOKKOS_INLINE_FUNCTION
void get_field(const int xz, const int isig_step,
const ScratchViewType<type_real> field,
const ScratchViewType<type_real> field_dot,
const ScratchViewType<type_real> field_dot_dot,
const ScratchViewType<type_real> hprime_xx,
const ScratchViewType<type_real> hprime_zz) const override;
void
get_field(const int xz, const int isig_step,
const ScratchViewType<type_real, medium::components> field,
const ScratchViewType<type_real, medium::components> field_dot,
const ScratchViewType<type_real, medium::components> field_dot_dot,
const ScratchViewType<type_real, 1> hprime_xx,
const ScratchViewType<type_real, 1> hprime_zz) const override;

KOKKOS_INLINE_FUNCTION
void compute_seismogram_components(
Expand Down
12 changes: 6 additions & 6 deletions include/domain/impl/receivers/acoustic/acoustic2d_isotropic.tpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,11 @@ KOKKOS_INLINE_FUNCTION void specfem::domain::impl::receivers::receiver<
specfem::enums::element::quadrature::static_quadrature_points<NGLL>,
specfem::enums::element::property::isotropic>::
get_field(const int xz, const int isig_step,
const ScratchViewType<type_real> field,
const ScratchViewType<type_real> field_dot,
const ScratchViewType<type_real> field_dot_dot,
const ScratchViewType<type_real> hprime_xx,
const ScratchViewType<type_real> hprime_zz) const {
const ScratchViewType<type_real, medium::components> field,
const ScratchViewType<type_real, medium::components> field_dot,
const ScratchViewType<type_real, medium::components> field_dot_dot,
const ScratchViewType<type_real, 1> hprime_xx,
const ScratchViewType<type_real, 1> hprime_zz) const {

#ifndef NDEBUG
assert(field.extent(0) == NGLL);
Expand All @@ -107,7 +107,7 @@ KOKKOS_INLINE_FUNCTION void specfem::domain::impl::receivers::receiver<
const type_real gammazl = this->gammaz(iz, ix);
const type_real rho_inversel = this->rho_inverse(iz, ix);

ScratchViewType<type_real> active_field;
ScratchViewType<type_real, 1> active_field;

switch (this->seismogram) {
case specfem::enums::seismogram::type::displacement:
Expand Down
17 changes: 10 additions & 7 deletions include/domain/impl/receivers/elastic/elastic2d.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,15 +17,18 @@ class receiver<specfem::enums::element::dimension::dim2,
using dimension = specfem::enums::element::dimension::dim2;
using medium = specfem::enums::element::medium::elastic;
using quadrature_points = qp_type;

template <typename T, int N>
using ScratchViewType =
typename quadrature_points::template ScratchViewType<type_real>;
typename quadrature_points::template ScratchViewType<T, N>;

KOKKOS_INLINE_FUNCTION virtual void get_field(
const int xz, const int isig_step, const ScratchViewType fieldx,
const ScratchViewType fieldz, const ScratchViewType fieldx_dot,
const ScratchViewType fieldz_dot, const ScratchViewType fieldx_dot_dot,
const ScratchViewType fieldz_dot_dot, const ScratchViewType s_hprime_xx,
const ScratchViewType s_hprime_zz) const {};
KOKKOS_INLINE_FUNCTION virtual void
get_field(const int xz, const int isig_step,
const ScratchViewType<type_real, medium::components> fieldx,
const ScratchViewType<type_real, medium::components> field_dot,
const ScratchViewType<type_real, medium::components> fieldx_dot_dot,
const ScratchViewType<type_real, 1> s_hprime_xx,
const ScratchViewType<type_real, 1> s_hprime_zz) const {};
KOKKOS_INLINE_FUNCTION virtual void compute_seismogram_components(
const int xz, const int isig_step,
dimension::array_type<type_real> &l_seismogram_components) const {};
Expand Down
17 changes: 7 additions & 10 deletions include/domain/impl/receivers/elastic/elastic2d_isotropic.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,9 +48,9 @@ class receiver<
using quadrature_points =
specfem::enums::element::quadrature::static_quadrature_points<NGLL>;

template <typename T>
template <typename T, int N>
using ScratchViewType =
typename quadrature_points::template ScratchViewType<T>;
typename quadrature_points::template ScratchViewType<T, N>;

KOKKOS_FUNCTION receiver() = default;

Expand All @@ -65,14 +65,11 @@ class receiver<

KOKKOS_INLINE_FUNCTION void
get_field(const int xz, const int isig_step,
const ScratchViewType<type_real> fieldx,
const ScratchViewType<type_real> fieldz,
const ScratchViewType<type_real> fieldx_dot,
const ScratchViewType<type_real> fieldz_dot,
const ScratchViewType<type_real> fieldx_dot_dot,
const ScratchViewType<type_real> fieldz_dot_dot,
const ScratchViewType<type_real> s_hprime_xx,
const ScratchViewType<type_real> s_hprime_zz) const override;
const ScratchViewType<type_real, medium::components> field,
const ScratchViewType<type_real, medium::components> field_dot,
const ScratchViewType<type_real, medium::components> field_dot_dot,
const ScratchViewType<type_real, 1> s_hprime_xx,
const ScratchViewType<type_real, 1> s_hprime_zz) const override;

KOKKOS_INLINE_FUNCTION void compute_seismogram_components(
const int xz, const int isig_step,
Expand Down
Loading

0 comments on commit d30b73a

Please sign in to comment.