From 42a681fecaab53993644228b21e8ee4584e4432f Mon Sep 17 00:00:00 2001 From: Rohit Kakodkar Date: Thu, 24 Aug 2023 16:35:08 -0400 Subject: [PATCH] Removed non-merged implementations --- include/domain/acoustic/acoustic_domain.hpp | 242 ----- include/domain/acoustic/acoustic_domain.tpp | 910 ------------------- include/domain/acoustic/interface.hpp | 7 - include/domain/elastic/elastic_domain.hpp | 243 ----- include/domain/elastic/elastic_domain.tpp | 924 -------------------- include/domain/elastic/interface.hpp | 7 - 6 files changed, 2333 deletions(-) delete mode 100644 include/domain/acoustic/acoustic_domain.hpp delete mode 100644 include/domain/acoustic/acoustic_domain.tpp delete mode 100644 include/domain/acoustic/interface.hpp delete mode 100644 include/domain/elastic/elastic_domain.hpp delete mode 100644 include/domain/elastic/elastic_domain.tpp delete mode 100644 include/domain/elastic/interface.hpp diff --git a/include/domain/acoustic/acoustic_domain.hpp b/include/domain/acoustic/acoustic_domain.hpp deleted file mode 100644 index 0aa66b82..00000000 --- a/include/domain/acoustic/acoustic_domain.hpp +++ /dev/null @@ -1,242 +0,0 @@ -#ifndef _ACOUSTIC_DOMAIN_HPP -#define _ACOUSTIC_DOMAIN_HPP - -#include "compute/interface.hpp" -#include "domain/domain.hpp" -#include "domain/impl/elements/interface.hpp" -#include "domain/impl/receivers/interface.hpp" -#include "domain/impl/sources/interface.hpp" -#include "quadrature/interface.hpp" -#include "specfem_enums.hpp" -#include "specfem_setup.hpp" -#include - -namespace specfem { -namespace domain { - -/** - * @brief Specialized domain class for acoustic media - * - * @tparam qp_type class used to define the quadrature points either at compile - * time or run time - */ -template -class domain { - -public: - using dimension = specfem::enums::element::dimension::dim2; - /** - * @brief Get a view of field stored on the device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_field() const { - return this->field; - } - /** - * @brief Get a view of field stored on the host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_field() const { - return this->h_field; - } - /** - * @brief Get a view of derivate of field stored on device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_field_dot() const { - return this->field_dot; - } - /** - * @brief Get a view of derivative of field stored on host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_field_dot() const { - return this->h_field_dot; - } - /** - * @brief Get a view of double derivative of field stored on device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_field_dot_dot() const { - return this->field_dot_dot; - } - /** - * @brief Get a view of double derivative of field stored on host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_field_dot_dot() const { - return this->h_field_dot_dot; - } - /** - * @brief Get a view of inverse of mass matrix stored on device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_rmass_inverse() const { - return this->rmass_inverse; - } - /** - * @brief Get a view of inverse of mass matrix stored on host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_rmass_inverse() const { - return this->h_rmass_inverse; - } - - /** - * @brief Construct a new domain object - * - * @param ndim Number of dimensions - * @param nglob Total number of distinct quadrature points - * @param quadrature_points quadrature points to define compile time - * quadrature or runtime quadrature - * @param compute Compute struct used to store global index of each quadrature - * point - * @param material_properties Struct used to store material properties at each - * quadrature point - * @param partial_derivatives Struct used to store partial derivatives of - * shape functions at each quadrature - * @param compute_sources Struct used to store pre-computed source arrays at - * each quadrature point and source time function - * @param receivers Struct used to store pre-computed receiver arrays at each - * quadrature point - * @param quadx quadrature class in x direction - * @param quadz quadrature class in z direction - */ - domain(const int ndim, const int nglob, const qp_type &quadrature_points, - specfem::compute::compute *compute, - specfem::compute::properties material_properties, - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::sources compute_sources, - specfem::compute::receivers receivers, - specfem::quadrature::quadrature *quadx, - specfem::quadrature::quadrature *quadz); - /** - * @brief Compute interaction of stiffness matrix on acceleration - * - */ - void compute_stiffness_interaction(); - - /** - * @brief Divide the acceleration by the mass matrix - * - */ - void divide_mass_matrix(); - /** - * @brief Compute interaction of sources on acceleration - * - * @param timeval - */ - void compute_source_interaction(const type_real timeval); - /** - * @brief Sync displacements views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_field(specfem::sync::kind kind); - /** - * @brief Sync velocity views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_field_dot(specfem::sync::kind kind); - /** - * @brief Sync acceleration views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_field_dot_dot(specfem::sync::kind kind); - /** - * @brief Sync inverse of mass matrix views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_rmass_inverse(specfem::sync::kind kind); - /** - * @brief Compute seismograms at for all receivers at isig_step - * - * @param seismogram_types DeviceView of types of seismograms to be - * calculated - * @param isig_step timestep for seismogram calculation - */ - void compute_seismogram(const int isig_step); - -private: - specfem::kokkos::DeviceView2d - field; ///< View of field on Device - specfem::kokkos::HostMirror2d - h_field; ///< View of field on host - specfem::kokkos::DeviceView2d - field_dot; ///< View of derivative of - ///< field on Device - specfem::kokkos::HostMirror2d - h_field_dot; ///< View of derivative - ///< of field on host - specfem::kokkos::DeviceView2d - field_dot_dot; ///< View of second - ///< derivative of - ///< field on Device - specfem::kokkos::HostMirror2d - h_field_dot_dot; ///< View of second - ///< derivative of - ///< field on host - specfem::kokkos::DeviceView2d - rmass_inverse; ///< View of inverse - ///< of mass matrix on - ///< device - specfem::kokkos::HostMirror2d - h_rmass_inverse; ///< View of inverse - ///< of mass matrix - ///< on host - specfem::compute::compute *compute; ///< Pointer to compute struct used to - ///< store spectral element numbering - ///< mapping (ibool) - quadrature::quadrature *quadx; ///< Pointer to quadrature object in - ///< x-dimension - quadrature::quadrature *quadz; ///< Pointer to quadrature object in - ///< z-dimension - int nelem_domain; ///< Total number of elements in this domain - specfem::kokkos::DeviceView1d > > - elements; ///< Container to store pointer to every element inside - ///< this domain - specfem::kokkos::DeviceView1d > > - sources; ///< Container to store pointer to every source inside - ///< this domain - specfem::kokkos::DeviceView1d > > - receivers; ///< Container to store pointer to every receiver inside - ///< this domain - - qp_type quadrature_points; ///< Quadrature points to define compile time - ///< quadrature or runtime quadrature -}; - -} // namespace domain -} // namespace specfem - -#endif diff --git a/include/domain/acoustic/acoustic_domain.tpp b/include/domain/acoustic/acoustic_domain.tpp deleted file mode 100644 index 6643429a..00000000 --- a/include/domain/acoustic/acoustic_domain.tpp +++ /dev/null @@ -1,910 +0,0 @@ -#ifndef _ACOUSTIC_DOMAIN_TPP -#define _ACOUSTIC_DOMAIN_TPP - -#include "compute/interface.hpp" -#include "constants.hpp" -#include "domain/acoustic/acoustic_domain.hpp" -#include "domain/domain.hpp" -#include "domain/impl/elements/interface.hpp" -#include "domain/impl/sources/interface.hpp" -#include "domain/impl/receivers/interface.hpp" -#include "globals.h" -#include "kokkos_abstractions.h" -#include "quadrature/interface.hpp" -#include "specfem_enums.hpp" -#include "specfem_setup.hpp" -#include -#include - -namespace acoustic_tmp { - -template -using element_container = - typename specfem::domain::impl::elements::container; - -template -using element_type = typename specfem::domain::impl::elements::element< - specfem::enums::element::dimension::dim2, - specfem::enums::element::medium::acoustic, qp_type, traits...>; - -template -using source_container = - typename specfem::domain::impl::sources::container; - -template -using source_type = typename specfem::domain::impl::sources::source< - specfem::enums::element::dimension::dim2, - specfem::enums::element::medium::acoustic, qp_type, traits...>; - -template -using receiver_container = - typename specfem::domain::impl::receivers::container; - -template -using receiver_type = typename specfem::domain::impl::receivers::receiver< - specfem::enums::element::dimension::dim2, - specfem::enums::element::medium::acoustic, qp_type, traits...>; - -template -void initialize_views( - const int &nglob, - specfem::kokkos::DeviceView2d field, - specfem::kokkos::DeviceView2d field_dot, - specfem::kokkos::DeviceView2d field_dot_dot, - specfem::kokkos::DeviceView2d - rmass_inverse) { - - constexpr int components = medium::components; - - Kokkos::parallel_for( - "specfem::Domain::acoustic::initiaze_views", - specfem::kokkos::DeviceMDrange<2, Kokkos::Iterate::Left>( - { 0, 0 }, { nglob, components }), - KOKKOS_LAMBDA(const int iglob, const int idim) { - field(iglob, idim) = 0; - field_dot(iglob, idim) = 0; - field_dot_dot(iglob, idim) = 0; - rmass_inverse(iglob, idim) = 0; - }); -} - -template -void initialize_rmass_inverse( - const specfem::kokkos::DeviceView3d ibool, - const specfem::kokkos::DeviceView1d - ispec_type, - const specfem::kokkos::DeviceView1d wxgll, - const specfem::kokkos::DeviceView1d wzgll, - const specfem::kokkos::DeviceView3d rho, - const specfem::kokkos::DeviceView3d kappa, - const specfem::kokkos::DeviceView3d jacobian, - specfem::kokkos::DeviceView2d - rmass_inverse) { - // Compute the mass matrix - - constexpr int components = medium::components; - constexpr auto value = medium::value; - - const int nspec = ibool.extent(0); - const int ngllz = ibool.extent(1); - const int ngllx = ibool.extent(2); - - const int nglob = rmass_inverse.extent(0); - - specfem::kokkos::DeviceScatterView2d results( - rmass_inverse); - - Kokkos::parallel_for( - "specfem::Domain::acoustic::compute_mass_matrix", - specfem::kokkos::DeviceMDrange<3, Kokkos::Iterate::Left>( - { 0, 0, 0 }, { nspec, ngllz, ngllx }), - KOKKOS_LAMBDA(const int ispec, const int iz, const int ix) { - int iglob = ibool(ispec, iz, ix); - type_real kappal = kappa(ispec, iz, ix); - auto access = results.access(); - if (ispec_type(ispec) == value) { - for (int icomponent = 0; icomponent < components; icomponent++) { - access(iglob, icomponent) += - wxgll(ix) * wzgll(iz) * jacobian(ispec, iz, ix) / kappal; - } - } - }); - - Kokkos::Experimental::contribute(rmass_inverse, results); - - // invert the mass matrix - Kokkos::parallel_for( - "specfem::Domain::acoustic::Invert_mass_matrix", - specfem::kokkos::DeviceRange(0, nglob), KOKKOS_LAMBDA(const int iglob) { - if (rmass_inverse(iglob, 0) > 0.0) { - for (int icomponent = 0; icomponent < components; icomponent++) { - type_real rmass_inverse_l = rmass_inverse(iglob, icomponent); - rmass_inverse(iglob, icomponent) = - 1.0 / rmass_inverse(iglob, icomponent); - } - } else { - for (int icomponent = 0; icomponent < components; icomponent++) { - rmass_inverse(iglob, icomponent) = 1.0; - } - } - }); -} - -template -void instantialize_element( - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::properties properties, - specfem::kokkos::DeviceView1d ispec_domain, - specfem::kokkos::DeviceView1d< - acoustic_tmp::element_container > > - elements) { - - auto h_elements = Kokkos::create_mirror_view(elements); - - Kokkos::parallel_for( - "specfem::domain::acoustic_isotropic::allocate_memory", - specfem::kokkos::HostRange(0, h_elements.extent(0)), - KOKKOS_LAMBDA(const int i) { - h_elements(i).element = (element_type *) - Kokkos::kokkos_malloc(sizeof( - element_type)); - }); - - Kokkos::deep_copy(elements, h_elements); - - Kokkos::parallel_for( - "specfem::domain::acoustic_isotropic::instantialize_element", - specfem::kokkos::DeviceRange(0, ispec_domain.extent(0)), - KOKKOS_LAMBDA(const int i) { - const int ispec = ispec_domain(i); - auto &element = elements(ispec).element; - new (element) - element_type( - ispec, partial_derivatives, properties); - }); - - return; -} - -template -void initialize_sources( - const specfem::compute::properties &properties, - const specfem::compute::sources &compute_sources, - specfem::kokkos::DeviceView1d > > - &sources) { - - const auto h_ispec_type = properties.h_ispec_type; - const auto ispec_array = compute_sources.ispec_array; - const auto h_ispec_array = compute_sources.h_ispec_array; - const auto kappa = properties.kappa; - - int nsources_domain = 0; - - // Count the number of sources in the domain - for (int isource = 0; isource < ispec_array.extent(0); isource++) { - if (h_ispec_type(h_ispec_array(isource)) == - specfem::enums::element::acoustic) { - nsources_domain++; - } - } - - sources = - specfem::kokkos::DeviceView1d > >( - "specfem::domain::acoustic_isotropic::sources", nsources_domain); - - specfem::kokkos::DeviceView1d my_sources( - "specfem::domain::acoustic_isotropic::my_sources", nsources_domain); - - auto h_my_sources = Kokkos::create_mirror_view(my_sources); - auto h_sources = Kokkos::create_mirror_view(sources); - - // Check if the source is in the domain - int index = 0; - for (int isource = 0; isource < ispec_array.extent(0); isource++) { - if (h_ispec_type(h_ispec_array(isource)) == - specfem::enums::element::acoustic) { - h_my_sources(index) = isource; - index++; - } - } - -#ifndef NDEBUG - assert(index == nsources_domain); -#endif - - Kokkos::deep_copy(my_sources, h_my_sources); - - // Allocate memory for the sources on the device - Kokkos::parallel_for( - "specfem::domain::acoustic_isotropic::allocate_memory", - specfem::kokkos::HostRange(0, nsources_domain), - KOKKOS_LAMBDA(const int isource) { - h_sources(isource).source = (source_type *) - Kokkos::kokkos_malloc(sizeof( - source_type)); - }); - - Kokkos::deep_copy(sources, h_sources); - - // Initialize the sources - Kokkos::parallel_for( - "specfem::domain::acoustic_isotropic::initialize_source", - specfem::kokkos::DeviceRange(0, nsources_domain), - KOKKOS_LAMBDA(const int isource) { - const int ispec = ispec_array(my_sources(isource)); - - specfem::forcing_function::stf *source_time_function = - compute_sources.stf_array(my_sources(isource)).T; - - const auto sv_kappa = - Kokkos::subview(kappa, ispec, Kokkos::ALL(), Kokkos::ALL()); - - const auto sv_source_array = - Kokkos::subview(compute_sources.source_array, my_sources(isource), - Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - - auto &source = sources(isource).source; - new (source) - source_type( - ispec, sv_kappa, sv_source_array, source_time_function); - }); - - Kokkos::fence(); - return; -} - -template -void initialize_receivers( - const specfem::compute::receivers &compute_receivers, - const specfem::compute::partial_derivatives &partial_derivatives, - const specfem::compute::properties &properties, - specfem::kokkos::DeviceView1d > > - &receivers) { - - const auto h_ispec_type = properties.h_ispec_type; - const auto ispec_array = compute_receivers.ispec_array; - const auto h_ispec_array = compute_receivers.h_ispec_array; - const auto seis_types = compute_receivers.seismogram_types; - - int nreceivers_domain = 0; - - // Count the number of receivers in the domain - for (int irec = 0; irec < ispec_array.extent(0); irec++) { - if (h_ispec_type(h_ispec_array(irec)) == - specfem::enums::element::acoustic) { - nreceivers_domain++; - } - } - - nreceivers_domain = nreceivers_domain * seis_types.extent(0); - - receivers = specfem::kokkos::DeviceView1d< - receiver_container > >( - "specfem::domain::acoustic_isotropic::receivers", nreceivers_domain); - - specfem::kokkos::DeviceView1d my_receivers( - "specfem::domain::acoustic_isotropic::my_receivers", nreceivers_domain); - - auto h_my_receivers = Kokkos::create_mirror_view(my_receivers); - auto h_receivers = Kokkos::create_mirror_view(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) { - h_my_receivers(index) = irec; - index++; - } - } - -#ifndef NDEBUG - assert(index * seis_types.extent(0) == nreceivers_domain); -#endif - - Kokkos::deep_copy(my_receivers, h_my_receivers); - - // Allocate memory for the receivers on the device - Kokkos::parallel_for( - "specfem::domain::acoustic_isotropic::allocate_memory", - specfem::kokkos::HostRange(0, nreceivers_domain), - KOKKOS_LAMBDA(const int irec) { - h_receivers(irec).receiver = (receiver_type *) - Kokkos::kokkos_malloc(sizeof( - receiver_type)); - }); - - Kokkos::deep_copy(receivers, h_receivers); - - Kokkos::parallel_for( - "specfem::domain::acoustic_isotropic::initialize_receiver", - specfem::kokkos::DeviceRange(0, nreceivers_domain), - KOKKOS_LAMBDA(const int inum) { - const int irec = my_receivers(inum / seis_types.extent(0)); - const int iseis = inum % seis_types.extent(0); - const int ispec = ispec_array(irec); - const auto seis_type = seis_types(iseis); - - const type_real cos_rec = compute_receivers.cos_recs(irec); - const type_real sin_rec = compute_receivers.sin_recs(irec); - - auto sv_receiver_array = - Kokkos::subview(compute_receivers.receiver_array, irec, Kokkos::ALL, - Kokkos::ALL, Kokkos::ALL); - - // Subview the seismogram array at current receiver and current - // seismogram value - auto sv_receiver_seismogram = - Kokkos::subview(compute_receivers.seismogram, Kokkos::ALL, iseis, - irec, Kokkos::ALL); - auto sv_receiver_field = - Kokkos::subview(compute_receivers.receiver_field, Kokkos::ALL, irec, - iseis, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL); - - auto &receiver = receivers(inum).receiver; - new (receiver) - receiver_type( - ispec, sin_rec, cos_rec, seis_type, sv_receiver_array, - sv_receiver_seismogram, partial_derivatives, properties, - sv_receiver_field); - }); - - return; -} - -template -void assign_elemental_properties( - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::properties properties, - specfem::kokkos::DeviceView2d field_dot_dot, - specfem::kokkos::DeviceView1d > > - &elements, - const int &nspec, const int &ngllz, const int &ngllx, int &nelem_domain) { - - // Assign elemental properties - // ---------------------------------------------------------------------- - nelem_domain = 0; - for (int ispec = 0; ispec < nspec; ispec++) { - if (properties.h_ispec_type(ispec) == - specfem::enums::element::medium::acoustic::value) { - nelem_domain++; - } - } - - specfem::kokkos::DeviceView1d ispec_domain( - "specfem::domain::acoustic::h_ispec_domain", nelem_domain); - specfem::kokkos::HostMirror1d h_ispec_domain = - Kokkos::create_mirror_view(ispec_domain); - - elements = specfem::kokkos::DeviceView1d< - acoustic_tmp::element_container > >( - "specfem::domain::acoustic::elements", nelem_domain); - - int index = 0; - for (int ispec = 0; ispec < nspec; ispec++) { - if (properties.h_ispec_type(ispec) == specfem::enums::element::acoustic) { - h_ispec_domain(index) = ispec; - index++; - } - } - - Kokkos::deep_copy(ispec_domain, h_ispec_domain); - - instantialize_element(partial_derivatives, properties, ispec_domain, - elements); -}; -} // namespace acoustic_tmp - -template -specfem::domain::domain:: - domain(const int ndim, const int nglob, const qp_type &quadrature_points, - specfem::compute::compute *compute, - specfem::compute::properties material_properties, - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::sources compute_sources, - specfem::compute::receivers compute_receivers, - specfem::quadrature::quadrature *quadx, - specfem::quadrature::quadrature *quadz) - : field(specfem::kokkos::DeviceView2d( - "specfem::Domain::acoustic::field", nglob, - specfem::enums::element::medium::acoustic::components)), - field_dot(specfem::kokkos::DeviceView2d( - "specfem::Domain::acoustic::field_dot", nglob, - specfem::enums::element::medium::acoustic::components)), - field_dot_dot( - specfem::kokkos::DeviceView2d( - "specfem::Domain::acoustic::field_dot_dot", nglob, - specfem::enums::element::medium::acoustic::components)), - rmass_inverse( - specfem::kokkos::DeviceView2d( - "specfem::Domain::acoustic::rmass_inverse", nglob, - specfem::enums::element::medium::acoustic::components)), - quadrature_points(quadrature_points), compute(compute), quadx(quadx), - quadz(quadz) { - - this->h_field = Kokkos::create_mirror_view(this->field); - this->h_field_dot = Kokkos::create_mirror_view(this->field_dot); - this->h_field_dot_dot = Kokkos::create_mirror_view(this->field_dot_dot); - this->h_rmass_inverse = Kokkos::create_mirror_view(this->rmass_inverse); - - const auto ibool = compute->ibool; - const int nspec = ibool.extent(0); - const int ngllz = ibool.extent(1); - const int ngllx = ibool.extent(2); - - //---------------------------------------------------------------------------- - // Initialize views - - acoustic_tmp::initialize_views( - nglob, this->field, this->field_dot, this->field_dot_dot, - this->rmass_inverse); - - //---------------------------------------------------------------------------- - // Inverse of mass matrix - - acoustic_tmp::initialize_rmass_inverse< - specfem::enums::element::medium::acoustic>( - compute->ibool, material_properties.ispec_type, quadx->get_w(), - quadz->get_w(), material_properties.rho, material_properties.kappa, - partial_derivatives.jacobian, this->rmass_inverse); - - // ---------------------------------------------------------------------------- - // Inverse of mass matrix - - acoustic_tmp::assign_elemental_properties( - partial_derivatives, material_properties, this->field_dot_dot, - this->elements, nspec, ngllz, ngllx, this->nelem_domain); - - // ---------------------------------------------------------------------------- - // Initialize the sources - - acoustic_tmp::initialize_sources(material_properties, compute_sources, this->sources); - - // ---------------------------------------------------------------------------- - // Initialize the receivers - acoustic_tmp::initialize_receivers(compute_receivers, partial_derivatives, - material_properties, this->receivers); - - return; -}; - -template -void specfem::domain::domain::sync_field(specfem::sync::kind kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_field, field); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(field, h_field); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::sync_field_dot(specfem::sync::kind - kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_field_dot, field_dot); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(field_dot, h_field_dot); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::sync_field_dot_dot(specfem::sync::kind - kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_field_dot_dot, field_dot_dot); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(field_dot_dot, h_field_dot_dot); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::sync_rmass_inverse(specfem::sync::kind - kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_rmass_inverse, rmass_inverse); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(rmass_inverse, h_rmass_inverse); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::divide_mass_matrix() { - - const int nglob = this->rmass_inverse.extent(0); - - Kokkos::parallel_for( - "specfem::Domain::acoustic::divide_mass_matrix", - specfem::kokkos::DeviceRange( - 0, specfem::enums::element::medium::acoustic::components * nglob), - KOKKOS_CLASS_LAMBDA(const int in) { - const int iglob = in % nglob; - const int idim = in / nglob; - this->field_dot_dot(iglob, idim) = - this->field_dot_dot(iglob, idim) * this->rmass_inverse(iglob, idim); - }); - - // Kokkos::fence(); - - return; -} - -template -void specfem::domain::domain::compute_stiffness_interaction() { - - const auto hprime_xx = this->quadx->get_hprime(); - const auto hprime_zz = this->quadz->get_hprime(); - const auto wxgll = this->quadx->get_w(); - const auto wzgll = this->quadz->get_w(); - const auto ibool = this->compute->ibool; - - int scratch_size = - 2 * - quadrature_points.template shmem_size(); - - scratch_size += - 2 * - quadrature_points.template shmem_size(); - - scratch_size += - 3 * - quadrature_points.template shmem_size(); - - scratch_size += - quadrature_points.template shmem_size(); - - Kokkos::parallel_for( - "specfem::Domain::acoustic::compute_gradients", - specfem::kokkos::DeviceTeam(this->nelem_domain, NTHREADS, NLANES) - .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 auto element = this->elements(team_member.league_rank()); - const auto ispec = element.get_ispec(); - - // Instantiate shared views - // --------------------------------------------------------------- - auto s_hprime_xx = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::x, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_hprime_zz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::z>( - team_member.team_scratch(0)); - auto s_hprimewgll_xx = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::x, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_hprimewgll_zz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::z>( - team_member.team_scratch(0)); - - auto s_field_chi = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_stress_integrand_xi = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_stress_integrand_gamma = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_iglob = - quadrature_points.template ScratchView( - team_member.team_scratch(0)); - - // ---------- Allocate shared views ------------------------------- - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - s_hprime_xx(iz, ix) = hprime_xx(iz, ix); - s_hprimewgll_xx(ix, iz) = wxgll(iz) * hprime_xx(iz, ix); - }); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllz, iz, ix); - s_hprime_zz(iz, ix) = hprime_zz(iz, ix); - s_hprimewgll_zz(ix, iz) = wzgll(iz) * hprime_zz(iz, ix); - }); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - const int iglob = ibool(ispec, iz, ix); - s_field_chi(iz, ix) = field(iglob, 0); - s_stress_integrand_xi(iz, ix) = 0.0; - s_stress_integrand_gamma(iz, ix) = 0.0; - s_iglob(iz, ix) = iglob; - }); - - // ------------------------------------------------------------------ - - team_member.team_barrier(); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int ix, iz; - sub2ind(xz, ngllx, iz, ix); - - type_real dchidxl = 0.0; - type_real dchidzl = 0.0; - - element.compute_gradient(xz, s_hprime_xx, s_hprime_zz, - s_field_chi, &dchidxl, &dchidzl); - - element.compute_stress(xz, dchidxl, dchidzl, - &s_stress_integrand_xi(iz, ix), - &s_stress_integrand_gamma(iz, ix)); - }); - - team_member.team_barrier(); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - - const int iglob = s_iglob(iz, ix); - const type_real wxglll = wxgll(ix); - const type_real wzglll = wzgll(iz); - - auto sv_field_dot_dot = - Kokkos::subview(field_dot_dot, iglob, Kokkos::ALL()); - - element.update_acceleration( - xz, wxglll, wzglll, s_stress_integrand_xi, - s_stress_integrand_gamma, s_hprimewgll_xx, s_hprimewgll_zz, - sv_field_dot_dot); - }); - }); - - Kokkos::fence(); - - return; -} - -template -void specfem::domain::domain< - specfem::enums::element::medium::acoustic, - qp_type>::compute_source_interaction(const type_real timeval) { - - const int nsources = this->sources.extent(0); - const auto ibool = this->compute->ibool; - - Kokkos::parallel_for( - "specfem::Domain::acoustic::compute_source_interaction", - specfem::kokkos::DeviceTeam(nsources, Kokkos::AUTO, 1), - KOKKOS_CLASS_LAMBDA( - const specfem::kokkos::DeviceTeam::member_type &team_member) { - int ngllx, ngllz; - quadrature_points.get_ngll(&ngllx, &ngllz); - int isource = team_member.league_rank(); - const auto source = this->sources(isource); - const int ispec = source.get_ispec(); - auto sv_ibool = Kokkos::subview(ibool, ispec, Kokkos::ALL, Kokkos::ALL); - - type_real stf; - - Kokkos::parallel_reduce( - Kokkos::TeamThreadRange(team_member, 1), - [=](const int &, type_real &lsum) { - lsum = source.eval_stf(timeval); - }, - stf); - - team_member.team_barrier(); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [=](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - int iglob = ibool(ispec, iz, ix); - - type_real accel; - auto sv_field_dot_dot = - Kokkos::subview(field_dot_dot, iglob, Kokkos::ALL()); - - source.compute_interaction(xz, stf, &accel); - source.update_acceleration(accel, sv_field_dot_dot); - }); - }); - - Kokkos::fence(); - return; -} - -// TODO : Implement acoustic element compute seismogram -template -void specfem::domain::domain::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. - - 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(); - - // hprime_zz - scratch_size += - quadrature_points.template shmem_size(); - - // field, field_dot, field_dot_dot - x and z components - scratch_size += - 3 * - quadrature_points.template shmem_size(); - - Kokkos::parallel_for( - "specfem::Domain::acoustic::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, specfem::enums::axes::x, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_hprime_zz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::z>( - team_member.team_scratch(0)); - - auto s_field = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_field_dot = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_field_dot_dot = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - // Allocate shared views - // ---------------------------------------------------------------- - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - 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( - 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( - team_member), - [=](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - int iglob = ibool(ispec, iz, ix); - s_field(iz, ix) = this->field(iglob, 0); - s_field_dot(iz, ix) = this->field_dot(iglob, 0); - s_field_dot_dot(iz, ix) = this->field_dot_dot(iglob, 0); - }); - - // Get seismogram field - // ---------------------------------------------------------------- - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - 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 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 &l_seismogram_components) { - receiver.compute_seismogram_components(xz, isig_step, - l_seismogram_components); - }, - specfem::kokkos::Sum >( - seismogram_components)); - Kokkos::single(Kokkos::PerTeam(team_member), [=] { - receiver.compute_seismogram(isig_step, seismogram_components); - }); - break; - } - }); -} - -#endif diff --git a/include/domain/acoustic/interface.hpp b/include/domain/acoustic/interface.hpp deleted file mode 100644 index c716579c..00000000 --- a/include/domain/acoustic/interface.hpp +++ /dev/null @@ -1,7 +0,0 @@ -#ifndef _ACOUSTIC_DOMAIN_INTERFACE_HPP -#define _ACOUSTIC_DOMAIN_INTERFACE_HPP - -#include "acoustic_domain.hpp" -#include "acoustic_domain.tpp" - -#endif diff --git a/include/domain/elastic/elastic_domain.hpp b/include/domain/elastic/elastic_domain.hpp deleted file mode 100644 index 115b2b07..00000000 --- a/include/domain/elastic/elastic_domain.hpp +++ /dev/null @@ -1,243 +0,0 @@ -#ifndef _ELASTIC_DOMAIN_HPP -#define _ELASTIC_DOMAIN_HPP - -#include "compute/interface.hpp" -#include "domain/domain.hpp" -#include "domain/impl/elements/interface.hpp" -#include "domain/impl/receivers/interface.hpp" -#include "domain/impl/sources/interface.hpp" -#include "quadrature/interface.hpp" -#include "specfem_enums.hpp" -#include "specfem_setup.hpp" -#include - -namespace specfem { -namespace domain { - -/** - * @brief Specialized domain class for elastic media - * - * @tparam qp_type class used to define the quadrature points either at compile - * time or run time - */ -template -class domain { - -public: - using dimension = specfem::enums::element::dimension::dim2; - /** - * @brief Get a view of field stored on the device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_field() const { - return this->field; - } - /** - * @brief Get a view of field stored on the host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_field() const { - return this->h_field; - } - /** - * @brief Get a view of derivate of field stored on device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_field_dot() const { - return this->field_dot; - } - /** - * @brief Get a view of derivative of field stored on host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_field_dot() const { - return this->h_field_dot; - } - /** - * @brief Get a view of double derivative of field stored on device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_field_dot_dot() const { - return this->field_dot_dot; - } - /** - * @brief Get a view of double derivative of field stored on host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_field_dot_dot() const { - return this->h_field_dot_dot; - } - /** - * @brief Get a view of inverse of mass matrix stored on device - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::DeviceView2d - get_rmass_inverse() const { - return this->rmass_inverse; - } - /** - * @brief Get a view of inverse of mass matrix stored on host - * - * @return specfem::kokkos::DeviceView2d - */ - specfem::kokkos::HostMirror2d - get_host_rmass_inverse() const { - return this->h_rmass_inverse; - } - - /** - * @brief Construct a new domain object - * - * @param ndim Number of dimensions - * @param nglob Total number of distinct quadrature points - * @param quadrature_points quadrature points to define compile time - * quadrature or runtime quadrature - * @param compute Compute struct used to store global index of each quadrature - * point - * @param material_properties Struct used to store material properties at each - * quadrature point - * @param partial_derivatives Struct used to store partial derivatives of - * shape functions at each quadrature - * @param compute_sources Struct used to store pre-computed source arrays at - * each quadrature point and source time function - * @param receivers Struct used to store pre-computed receiver arrays at each - * quadrature point - * @param quadx quadrature class in x direction - * @param quadz quadrature class in z direction - */ - domain(const int ndim, const int nglob, const qp_type &quadrature_points, - specfem::compute::compute *compute, - specfem::compute::properties material_properties, - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::sources compute_sources, - specfem::compute::receivers compute_receivers, - specfem::quadrature::quadrature *quadx, - specfem::quadrature::quadrature *quadz); - /** - * @brief Compute interaction of stiffness matrix on acceleration - * - */ - void compute_stiffness_interaction(); - - /** - * @brief Divide the acceleration by the mass matrix - * - */ - void divide_mass_matrix(); - /** - * @brief Compute interaction of sources on acceleration - * - * @param timeval - */ - void compute_source_interaction(const type_real timeval); - /** - * @brief Sync displacements views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_field(specfem::sync::kind kind); - /** - * @brief Sync velocity views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_field_dot(specfem::sync::kind kind); - /** - * @brief Sync acceleration views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_field_dot_dot(specfem::sync::kind kind); - /** - * @brief Sync inverse of mass matrix views between host and device - * - * @param kind defines sync direction i.e. DeviceToHost or HostToDevice - */ - void sync_rmass_inverse(specfem::sync::kind kind); - /** - * @brief Compute seismograms at for all receivers at isig_step - * - * @param seismogram_types DeviceView of types of seismograms to be - * calculated - * @param isig_step timestep for seismogram calculation - */ - void compute_seismogram(const int isig_step); - -private: - specfem::kokkos::DeviceView2d - field; ///< View of field on Device - specfem::kokkos::HostMirror2d - h_field; ///< View of field on host - specfem::kokkos::DeviceView2d - field_dot; ///< View of derivative of - ///< field on Device - specfem::kokkos::HostMirror2d - h_field_dot; ///< View of derivative - ///< of field on host - specfem::kokkos::DeviceView2d - field_dot_dot; ///< View of second - ///< derivative of - ///< field on Device - specfem::kokkos::HostMirror2d - h_field_dot_dot; ///< View of second - ///< derivative of - ///< field on host - specfem::kokkos::DeviceView2d - rmass_inverse; ///< View of inverse - ///< of mass matrix on - ///< device - specfem::kokkos::HostMirror2d - h_rmass_inverse; ///< View of inverse - ///< of mass matrix - ///< on host - specfem::compute::compute *compute; ///< Pointer to compute struct used to - ///< store spectral element numbering - ///< mapping (ibool) - quadrature::quadrature *quadx; ///< Pointer to quadrature object in - ///< x-dimension - quadrature::quadrature *quadz; ///< Pointer to quadrature object in - ///< z-dimension - int nelem_domain; ///< Total number of elements in this domain - specfem::kokkos::DeviceView1d > > - elements; ///< Container to store pointer to every element inside - ///< this domain - specfem::kokkos::DeviceView1d > > - sources; ///< Container to store pointer to every source inside - ///< this domain - specfem::kokkos::DeviceView1d > > - receivers; ///< Container to store pointer to every receiver inside - ///< this domain - - qp_type quadrature_points; ///< Quadrature points to define compile time - ///< quadrature or runtime quadrature -}; - -} // namespace domain -} // namespace specfem - -#endif diff --git a/include/domain/elastic/elastic_domain.tpp b/include/domain/elastic/elastic_domain.tpp deleted file mode 100644 index 8b143fd1..00000000 --- a/include/domain/elastic/elastic_domain.tpp +++ /dev/null @@ -1,924 +0,0 @@ -#ifndef _ELASTIC_DOMAIN_TPP -#define _ELASTIC_DOMAIN_TPP - -#include "compute/interface.hpp" -#include "constants.hpp" -#include "domain/domain.hpp" -#include "domain/elastic/elastic_domain.hpp" -#include "domain/impl/elements/interface.hpp" -#include "domain/impl/receivers/interface.hpp" -#include "domain/impl/sources/interface.hpp" -#include "globals.h" -#include "kokkos_abstractions.h" -#include "quadrature/interface.hpp" -#include "specfem_enums.hpp" -#include "specfem_setup.hpp" -#include -#include - -// Type aliases -// ---------------------------------------------------------------------------- - -template -using element_container = - typename specfem::domain::impl::elements::container; - -template -using element_type = typename specfem::domain::impl::elements::element< - specfem::enums::element::dimension::dim2, - specfem::enums::element::medium::elastic, qp_type, traits...>; - -template -using source_container = - typename specfem::domain::impl::sources::container; - -template -using source_type = typename specfem::domain::impl::sources::source< - specfem::enums::element::dimension::dim2, - specfem::enums::element::medium::elastic, qp_type, traits...>; - -template -using receiver_container = - typename specfem::domain::impl::receivers::container; - -template -using receiver_type = typename specfem::domain::impl::receivers::receiver< - specfem::enums::element::dimension::dim2, - specfem::enums::element::medium::elastic, qp_type, traits...>; - -// ---------------------------------------------------------------------------- - -template -void initialize_views( - const int &nglob, - specfem::kokkos::DeviceView2d field, - specfem::kokkos::DeviceView2d field_dot, - specfem::kokkos::DeviceView2d field_dot_dot, - specfem::kokkos::DeviceView2d - rmass_inverse) { - - constexpr int components = medium::components; - - Kokkos::parallel_for( - "specfem::Domain::Elastic::initiaze_views", - specfem::kokkos::DeviceMDrange<2, Kokkos::Iterate::Left>( - { 0, 0 }, { nglob, components }), - KOKKOS_LAMBDA(const int iglob, const int idim) { - field(iglob, idim) = 0; - field_dot(iglob, idim) = 0; - field_dot_dot(iglob, idim) = 0; - rmass_inverse(iglob, idim) = 0; - }); -} - -template -void initialize_rmass_inverse( - const specfem::kokkos::DeviceView3d ibool, - const specfem::kokkos::DeviceView1d - ispec_type, - const specfem::kokkos::DeviceView1d wxgll, - const specfem::kokkos::DeviceView1d wzgll, - const specfem::kokkos::DeviceView3d rho, - const specfem::kokkos::DeviceView3d jacobian, - specfem::kokkos::DeviceView2d - rmass_inverse) { - // Compute the mass matrix - - constexpr int components = medium::components; - constexpr auto value = medium::value; - - const int nspec = ibool.extent(0); - const int ngllz = ibool.extent(1); - const int ngllx = ibool.extent(2); - - const int nglob = rmass_inverse.extent(0); - - specfem::kokkos::DeviceScatterView2d results( - rmass_inverse); - - Kokkos::parallel_for( - "specfem::Domain::Elastic::compute_mass_matrix", - specfem::kokkos::DeviceMDrange<3, Kokkos::Iterate::Left>( - { 0, 0, 0 }, { nspec, ngllz, ngllx }), - KOKKOS_LAMBDA(const int ispec, const int iz, const int ix) { - int iglob = ibool(ispec, iz, ix); - type_real rhol = rho(ispec, iz, ix); - auto access = results.access(); - if (ispec_type(ispec) == value) { - for (int icomponent = 0; icomponent < components; icomponent++) { - access(iglob, icomponent) += - wxgll(ix) * wzgll(iz) * rhol * jacobian(ispec, iz, ix); - } - } - }); - - Kokkos::Experimental::contribute(rmass_inverse, results); - - // invert the mass matrix - Kokkos::parallel_for( - "specfem::Domain::Elastic::Invert_mass_matrix", - specfem::kokkos::DeviceRange(0, nglob), KOKKOS_LAMBDA(const int iglob) { - if (rmass_inverse(iglob, 0) > 0.0) { - for (int icomponent = 0; icomponent < components; icomponent++) { - rmass_inverse(iglob, icomponent) = - 1.0 / rmass_inverse(iglob, icomponent); - } - } else { - for (int icomponent = 0; icomponent < components; icomponent++) { - rmass_inverse(iglob, icomponent) = 1.0; - } - } - }); -} - -template -void instantialize_element( - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::properties properties, - specfem::kokkos::DeviceView1d ispec_domain, - specfem::kokkos::DeviceView1d > > - elements) { - auto h_elements = Kokkos::create_mirror_view(elements); - - Kokkos::parallel_for( - "specfem::domain::elastic_isotropic::allocate_memory", - specfem::kokkos::HostRange(0, h_elements.extent(0)), - KOKKOS_LAMBDA(const int i) { - h_elements(i).element = (element_type *) - Kokkos::kokkos_malloc(sizeof( - element_type)); - }); - - Kokkos::deep_copy(elements, h_elements); - - Kokkos::parallel_for( - "specfem::domain::elastic_isotropic::instantialize_element", - specfem::kokkos::DeviceRange(0, ispec_domain.extent(0)), - KOKKOS_LAMBDA(const int i) { - const int ispec = ispec_domain(i); - auto &element = elements(ispec).element; - new (element) - element_type( - ispec, partial_derivatives, properties); - }); - - return; -} - -template -void initialize_sources( - const specfem::compute::properties &properties, - const specfem::compute::sources &compute_sources, - specfem::kokkos::DeviceView1d > > - &sources) { - - const auto h_ispec_type = properties.h_ispec_type; - const auto ispec_type = properties.ispec_type; - const auto h_ispec_array = compute_sources.h_ispec_array; - const auto ispec_array = compute_sources.ispec_array; - int nsources_domain = 0; - - specfem::kokkos::DeviceView1d my_recs; - - // Check how many sources belong to this domain - for (int ispec = 0; ispec < h_ispec_array.extent(0); ispec++) { - if (h_ispec_type(h_ispec_array(ispec)) == - specfem::enums::element::elastic) { - nsources_domain++; - } - } - - sources = - specfem::kokkos::DeviceView1d > >( - "specfem::domain::elastic_isotropic::sources", nsources_domain); - - specfem::kokkos::DeviceView1d my_sources( - "specfem::domain::elastic_isotropic::my_sources", nsources_domain); - - auto h_sources = Kokkos::create_mirror_view(sources); - auto h_my_sources = Kokkos::create_mirror_view(my_sources); - - // Store which sources belong to this domain - int index = 0; - for (int isource = 0; isource < h_ispec_array.extent(0); isource++) { - if (h_ispec_type(h_ispec_array(isource)) == - specfem::enums::element::elastic) { - h_my_sources(index) = isource; - index++; - } - } - -#ifndef NDEBUG - assert(index == nsources_domain); -#endif - - Kokkos::deep_copy(my_sources, h_my_sources); - - // Allocate memory for the sources on the device - Kokkos::parallel_for( - "specfem::domain::elastic_isotropic::allocate_memory", - specfem::kokkos::HostRange(0, nsources_domain), - KOKKOS_LAMBDA(const int i) { - h_sources(i).source = (source_type *) - Kokkos::kokkos_malloc(sizeof( - source_type)); - }); - - Kokkos::deep_copy(sources, h_sources); - - Kokkos::parallel_for( - "specfem::domain::elastic_isotropic::initialize_source", - specfem::kokkos::DeviceRange(0, nsources_domain), - KOKKOS_LAMBDA(const int isource) { - auto &source = sources(isource).source; - const int ispec = ispec_array(my_sources(isource)); - auto source_time_function = - compute_sources.stf_array(my_sources(isource)).T; - const auto sv_source_array = - Kokkos::subview(compute_sources.source_array, my_sources(isource), - Kokkos::ALL(), Kokkos::ALL(), Kokkos::ALL()); - new (source) - source_type( - ispec, sv_source_array, source_time_function); - }); - - return; -} - -template -void initialize_receivers( - const specfem::compute::receivers compute_receivers, - const specfem::compute::partial_derivatives &partial_derivatives, - const specfem::compute::properties &properties, - specfem::kokkos::DeviceView1d > > - &receivers) { - - const auto h_ispec_type = properties.h_ispec_type; - const auto h_ispec_array = compute_receivers.h_ispec_array; - const auto ispec_array = compute_receivers.ispec_array; - const auto seis_types = compute_receivers.seismogram_types; - - int nreceivers_domain = 0; - - // Check how many receivers belong to this domain - for (int irec = 0; irec < ispec_array.extent(0); irec++) { - if (h_ispec_type(h_ispec_array(irec)) == specfem::enums::element::elastic) { - nreceivers_domain++; - } - } - - nreceivers_domain = nreceivers_domain * seis_types.extent(0); - - receivers = specfem::kokkos::DeviceView1d< - receiver_container > >( - "specfem::domain::elastic_isotropic::receivers", nreceivers_domain); - - specfem::kokkos::DeviceView1d my_receivers( - "specfem::domain::elastic_isotropic::my_recs", nreceivers_domain); - - auto h_my_receivers = Kokkos::create_mirror_view(my_receivers); - auto h_receivers = Kokkos::create_mirror_view(receivers); - - // Store which receivers belong to this 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::elastic) { - h_my_receivers(index) = irec; - index++; - } - } - -#ifndef NDEBUG - assert(index * seis_types.extent(0) == nreceivers_domain); -#endif - - Kokkos::parallel_for( - "specfem::domain::elastic_isotropic::allocate_memory", - specfem::kokkos::HostRange(0, nreceivers_domain), - KOKKOS_LAMBDA(const int i) { - h_receivers(i).receiver = (receiver_type *) - Kokkos::kokkos_malloc(sizeof( - receiver_type)); - }); - - Kokkos::deep_copy(receivers, h_receivers); - Kokkos::deep_copy(my_receivers, h_my_receivers); - - Kokkos::parallel_for( - "specfem::domain::elastic_isotropic::initialize_receiver", - specfem::kokkos::DeviceRange(0, nreceivers_domain), - KOKKOS_LAMBDA(const int inum) { - const int irec = my_receivers(inum / seis_types.extent(0)); - const int iseis = inum % seis_types.extent(0); - - const int ispec = ispec_array(irec); - const auto seis_type = seis_types(iseis); - - const type_real cos_rec = compute_receivers.cos_recs(irec); - const type_real sin_rec = compute_receivers.sin_recs(irec); - - auto sv_receiver_array = - Kokkos::subview(compute_receivers.receiver_array, irec, Kokkos::ALL, - Kokkos::ALL, Kokkos::ALL); - - // Subview the seismogram array at current receiver and current - // seismogram value - auto sv_receiver_seismogram = - Kokkos::subview(compute_receivers.seismogram, Kokkos::ALL, iseis, - irec, Kokkos::ALL); - auto sv_receiver_field = - Kokkos::subview(compute_receivers.receiver_field, Kokkos::ALL, irec, - iseis, Kokkos::ALL, Kokkos::ALL, Kokkos::ALL); - - auto &receiver = receivers(inum).receiver; - new (receiver) - receiver_type( - ispec, sin_rec, cos_rec, seis_type, sv_receiver_array, - sv_receiver_seismogram, partial_derivatives, properties, - sv_receiver_field); - }); - - return; -} - -template -void assign_elemental_properties( - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::properties properties, - specfem::kokkos::DeviceView2d field_dot_dot, - specfem::kokkos::DeviceView1d > > - &elements, - const int &nspec, const int &ngllz, const int &ngllx, int &nelem_domain) { - - // Assign elemental properties - // ---------------------------------------------------------------------- - nelem_domain = 0; - for (int ispec = 0; ispec < nspec; ispec++) { - if (properties.h_ispec_type(ispec) == - specfem::enums::element::medium::elastic::value) { - nelem_domain++; - } - } - - specfem::kokkos::DeviceView1d ispec_domain( - "specfem::domain::elastic::h_ispec_domain", nelem_domain); - specfem::kokkos::HostMirror1d h_ispec_domain = - Kokkos::create_mirror_view(ispec_domain); - - elements = - specfem::kokkos::DeviceView1d > >( - "specfem::domain::elastic::elements", nelem_domain); - - int index = 0; - for (int ispec = 0; ispec < nspec; ispec++) { - if (properties.h_ispec_type(ispec) == specfem::enums::element::elastic) { - h_ispec_domain(index) = ispec; - index++; - } - } - - Kokkos::deep_copy(ispec_domain, h_ispec_domain); - - instantialize_element(partial_derivatives, properties, ispec_domain, - elements); -}; - -template -specfem::domain::domain:: - domain(const int ndim, const int nglob, const qp_type &quadrature_points, - specfem::compute::compute *compute, - specfem::compute::properties material_properties, - specfem::compute::partial_derivatives partial_derivatives, - specfem::compute::sources compute_sources, - specfem::compute::receivers compute_receivers, - specfem::quadrature::quadrature *quadx, - specfem::quadrature::quadrature *quadz) - : field(specfem::kokkos::DeviceView2d( - "specfem::Domain::Elastic::field", nglob, - specfem::enums::element::medium::elastic::components)), - field_dot(specfem::kokkos::DeviceView2d( - "specfem::Domain::Elastic::field_dot", nglob, - specfem::enums::element::medium::elastic::components)), - field_dot_dot( - specfem::kokkos::DeviceView2d( - "specfem::Domain::Elastic::field_dot_dot", nglob, - specfem::enums::element::medium::elastic::components)), - rmass_inverse( - specfem::kokkos::DeviceView2d( - "specfem::Domain::Elastic::rmass_inverse", nglob, - specfem::enums::element::medium::elastic::components)), - quadrature_points(quadrature_points), compute(compute), quadx(quadx), - quadz(quadz) { - - this->h_field = Kokkos::create_mirror_view(this->field); - this->h_field_dot = Kokkos::create_mirror_view(this->field_dot); - this->h_field_dot_dot = Kokkos::create_mirror_view(this->field_dot_dot); - this->h_rmass_inverse = Kokkos::create_mirror_view(this->rmass_inverse); - - const auto ibool = compute->ibool; - const int nspec = ibool.extent(0); - const int ngllz = ibool.extent(1); - const int ngllx = ibool.extent(2); - - //---------------------------------------------------------------------------- - // Initialize views - - initialize_views( - nglob, this->field, this->field_dot, this->field_dot_dot, - this->rmass_inverse); - - //---------------------------------------------------------------------------- - // Inverse of mass matrix - - initialize_rmass_inverse( - compute->ibool, material_properties.ispec_type, quadx->get_w(), - quadz->get_w(), material_properties.rho, partial_derivatives.jacobian, - this->rmass_inverse); - - // ---------------------------------------------------------------------------- - // Inverse of mass matrix - - assign_elemental_properties(partial_derivatives, material_properties, - this->field_dot_dot, this->elements, nspec, ngllz, - ngllx, this->nelem_domain); - - // ---------------------------------------------------------------------------- - // Initialize the sources - - initialize_sources(material_properties, compute_sources, this->sources); - - // ---------------------------------------------------------------------------- - // Initialize the receivers - - initialize_receivers(compute_receivers, partial_derivatives, - material_properties, this->receivers); - - return; -}; - -template -void specfem::domain::domain::sync_field(specfem::sync::kind kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_field, field); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(field, h_field); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::sync_field_dot(specfem::sync::kind - kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_field_dot, field_dot); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(field_dot, h_field_dot); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::sync_field_dot_dot(specfem::sync::kind - kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_field_dot_dot, field_dot_dot); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(field_dot_dot, h_field_dot_dot); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::sync_rmass_inverse(specfem::sync::kind - kind) { - - if (kind == specfem::sync::DeviceToHost) { - Kokkos::deep_copy(h_rmass_inverse, rmass_inverse); - } else if (kind == specfem::sync::HostToDevice) { - Kokkos::deep_copy(rmass_inverse, h_rmass_inverse); - } else { - throw std::runtime_error("Could not recognize the kind argument"); - } - - return; -} - -template -void specfem::domain::domain::divide_mass_matrix() { - - const int nglob = this->rmass_inverse.extent(0); - - Kokkos::parallel_for( - "specfem::Domain::Elastic::divide_mass_matrix", - specfem::kokkos::DeviceRange( - 0, specfem::enums::element::medium::elastic::components * nglob), - KOKKOS_CLASS_LAMBDA(const int in) { - const int iglob = in % nglob; - const int idim = in / nglob; - this->field_dot_dot(iglob, idim) = - this->field_dot_dot(iglob, idim) * this->rmass_inverse(iglob, idim); - }); - - // Kokkos::fence(); - - return; -} - -template -void specfem::domain::domain::compute_stiffness_interaction() { - - const auto hprime_xx = this->quadx->get_hprime(); - const auto hprime_zz = this->quadz->get_hprime(); - const auto wxgll = this->quadx->get_w(); - const auto wzgll = this->quadz->get_w(); - const auto ibool = this->compute->ibool; - - int scratch_size = - 2 * - quadrature_points.template shmem_size(); - - scratch_size += - 2 * - quadrature_points.template shmem_size(); - - scratch_size += - 6 * - quadrature_points.template shmem_size(); - - scratch_size += - quadrature_points.template shmem_size(); - - Kokkos::parallel_for( - "specfem::Domain::Elastic::compute_gradients", - specfem::kokkos::DeviceTeam(this->nelem_domain, NTHREADS, NLANES) - .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 auto element = this->elements(team_member.league_rank()); - const auto ispec = element.get_ispec(); - - // Instantiate shared views - // --------------------------------------------------------------- - auto s_hprime_xx = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::x, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_hprime_zz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::z>( - team_member.team_scratch(0)); - auto s_hprimewgll_xx = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::x, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_hprimewgll_zz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::z>( - team_member.team_scratch(0)); - - auto s_fieldx = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_fieldz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_stress_integral_1 = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_stress_integral_2 = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_stress_integral_3 = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_stress_integral_4 = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - auto s_iglob = - quadrature_points.template ScratchView( - team_member.team_scratch(0)); - - // ---------- Allocate shared views ------------------------------- - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - s_hprime_xx(iz, ix) = hprime_xx(iz, ix); - s_hprimewgll_xx(ix, iz) = wxgll(iz) * hprime_xx(iz, ix); - }); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllz, iz, ix); - s_hprime_zz(iz, ix) = hprime_zz(iz, ix); - s_hprimewgll_zz(ix, iz) = wzgll(iz) * hprime_zz(iz, ix); - }); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - const int iglob = ibool(ispec, iz, ix); - s_fieldx(iz, ix) = field(iglob, 0); - s_fieldz(iz, ix) = field(iglob, 1); - s_stress_integral_1(iz, ix) = 0.0; - s_stress_integral_2(iz, ix) = 0.0; - s_stress_integral_3(iz, ix) = 0.0; - s_stress_integral_4(iz, ix) = 0.0; - s_iglob(iz, ix) = iglob; - }); - - // ------------------------------------------------------------------ - - team_member.team_barrier(); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int ix, iz; - sub2ind(xz, ngllx, iz, ix); - - type_real duxdxl = 0.0; - type_real duxdzl = 0.0; - type_real duzdxl = 0.0; - type_real duzdzl = 0.0; - - element.compute_gradient(xz, s_hprime_xx, s_hprime_zz, s_fieldx, - s_fieldz, &duxdxl, &duxdzl, &duzdxl, - &duzdzl); - - element.compute_stress( - xz, duxdxl, duxdzl, duzdxl, duzdzl, - &s_stress_integral_1(iz, ix), &s_stress_integral_2(iz, ix), - &s_stress_integral_3(iz, ix), &s_stress_integral_4(iz, ix)); - }); - - team_member.team_barrier(); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [&](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - - const int iglob = s_iglob(iz, ix); - const type_real wxglll = wxgll(ix); - const type_real wzglll = wzgll(iz); - - auto sv_field_dot_dot = - Kokkos::subview(field_dot_dot, iglob, Kokkos::ALL()); - - element.update_acceleration( - xz, wxglll, wzglll, s_stress_integral_1, s_stress_integral_2, - s_stress_integral_3, s_stress_integral_4, s_hprimewgll_xx, - s_hprimewgll_zz, sv_field_dot_dot); - }); - }); - return; -} - -template -void specfem::domain::domain< - specfem::enums::element::medium::elastic, - qp_type>::compute_source_interaction(const type_real timeval) { - - const int nsources = this->sources.extent(0); - const auto ibool = this->compute->ibool; - - Kokkos::parallel_for( - "specfem::Domain::Elastic::compute_source_interaction", - specfem::kokkos::DeviceTeam(nsources, Kokkos::AUTO, 1), - KOKKOS_CLASS_LAMBDA( - const specfem::kokkos::DeviceTeam::member_type &team_member) { - int ngllx, ngllz; - quadrature_points.get_ngll(&ngllx, &ngllz); - int isource = team_member.league_rank(); - const auto source = this->sources(isource); - const int ispec = source.get_ispec(); - auto sv_ibool = Kokkos::subview(ibool, ispec, Kokkos::ALL, Kokkos::ALL); - - type_real stf; - - Kokkos::parallel_reduce( - Kokkos::TeamThreadRange(team_member, 1), - [=](const int &, type_real &lsum) { - lsum = source.eval_stf(timeval); - }, - stf); - - team_member.team_barrier(); - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [=](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - int iglob = ibool(ispec, iz, ix); - - type_real accelx, accelz; - auto sv_field_dot_dot = - Kokkos::subview(field_dot_dot, iglob, Kokkos::ALL()); - - source.compute_interaction(xz, stf, &accelx, &accelz); - source.update_acceleration(accelx, accelz, sv_field_dot_dot); - }); - }); - return; -} - -template -void specfem::domain::domain::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. - - 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(); - - // hprime_zz - scratch_size += - quadrature_points.template shmem_size(); - - // field, field_dot, field_dot_dot - x and z components - scratch_size += - 6 * - quadrature_points.template shmem_size(); - - Kokkos::parallel_for( - "specfem::Domain::Elastic::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, specfem::enums::axes::x, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_hprime_zz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::z>( - team_member.team_scratch(0)); - - auto s_fieldx = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_fieldx_dot = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_fieldx_dot_dot = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_fieldz = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_fieldz_dot = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - auto s_fieldz_dot_dot = quadrature_points.template ScratchView< - type_real, specfem::enums::axes::z, specfem::enums::axes::x>( - team_member.team_scratch(0)); - - // Allocate shared views - // ---------------------------------------------------------------- - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - 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( - 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( - team_member), - [=](const int xz) { - int iz, ix; - sub2ind(xz, ngllx, iz, ix); - int iglob = ibool(ispec, iz, ix); - s_fieldx(iz, ix) = this->field(iglob, 0); - s_fieldx_dot(iz, ix) = this->field_dot(iglob, 0); - s_fieldx_dot_dot(iz, ix) = this->field_dot_dot(iglob, 0); - s_fieldz(iz, ix) = this->field(iglob, 1); - s_fieldz_dot(iz, ix) = this->field_dot(iglob, 1); - s_fieldz_dot_dot(iz, ix) = this->field_dot_dot(iglob, 1); - }); - - // Get seismogram field - // ---------------------------------------------------------------- - - Kokkos::parallel_for( - quadrature_points.template TeamThreadRange( - team_member), - [=](const int xz) { - receiver.get_field(xz, isig_step, s_fieldx, s_fieldz, - s_fieldx_dot, s_fieldz_dot, s_fieldx_dot_dot, - s_fieldz_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 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 &l_seismogram_components) { - receiver.compute_seismogram_components(xz, isig_step, - l_seismogram_components); - }, - specfem::kokkos::Sum >( - seismogram_components)); - Kokkos::single(Kokkos::PerTeam(team_member), [=] { - receiver.compute_seismogram(isig_step, seismogram_components); - }); - break; - } - }); -} - -#endif diff --git a/include/domain/elastic/interface.hpp b/include/domain/elastic/interface.hpp deleted file mode 100644 index cad736ea..00000000 --- a/include/domain/elastic/interface.hpp +++ /dev/null @@ -1,7 +0,0 @@ -#ifndef _ELASTIC_DOMAIN_INTERFACE_HPP -#define _ELASTIC_DOMAIN_INTERFACE_HPP - -#include "elastic_domain.hpp" -#include "elastic_domain.tpp" - -#endif