From 90129c964678d51c8aafe3e1a666d9946c3bbd1c Mon Sep 17 00:00:00 2001 From: Rohit Kakodkar Date: Wed, 27 Mar 2024 05:23:45 -0400 Subject: [PATCH] Added adjoint sources - Methods to read adjoint source file - User Defined sources --- CMakeLists.txt | 1 + .../compute/sources/impl/source_medium.tpp | 4 +- include/domain/impl/sources/kernel.tpp | 2 +- include/kokkos_abstractions.h | 3 + include/source/adjoint_source.hpp | 2 + include/source/external.hpp | 39 +++++++++ include/source/interface.hpp | 1 + include/source_time_function/external.hpp | 2 + src/reader/seismogram.cpp | 3 +- src/source/adjoint_source.cpp | 40 +++++++-- src/source/external.cpp | 83 +++++++++++++++++++ src/source/force_source.cpp | 2 +- src/source/moment_tensor_source.cpp | 2 +- src/source/read_sources.cpp | 6 ++ src/source_time_function/external.cpp | 19 +++-- 15 files changed, 189 insertions(+), 20 deletions(-) create mode 100644 include/source/external.hpp create mode 100644 src/source/external.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2b7a3f13..a1f23640 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -233,6 +233,7 @@ add_library( src/source/force_source.cpp src/source/moment_tensor_source.cpp src/source/adjoint_source.cpp + src/source/external.cpp src/source/read_sources.cpp ) diff --git a/include/compute/sources/impl/source_medium.tpp b/include/compute/sources/impl/source_medium.tpp index 7cfb348b..1e6271a5 100644 --- a/include/compute/sources/impl/source_medium.tpp +++ b/include/compute/sources/impl/source_medium.tpp @@ -22,8 +22,8 @@ specfem::compute::impl::sources::source_medium:: sources.size(), medium_type::components), h_source_time_function(Kokkos::create_mirror_view(source_time_function)), source_array("specfem::sources::source_array", sources.size(), - mesh.quadratures.gll.N, mesh.quadratures.gll.N, - medium_type::components), + medium_type::components, mesh.quadratures.gll.N, + mesh.quadratures.gll.N), h_source_array(Kokkos::create_mirror_view(source_array)) { for (int isource = 0; isource < sources.size(); isource++) { diff --git a/include/domain/impl/sources/kernel.tpp b/include/domain/impl/sources/kernel.tpp index 0d0a6ab0..a50e37af 100644 --- a/include/domain/impl/sources/kernel.tpp +++ b/include/domain/impl/sources/kernel.tpp @@ -98,7 +98,7 @@ void specfem::domain::impl::kernels::source_kernel< const specfem::kokkos::array_type lagrange_interpolant(Kokkos::subview( - sources.source_array, isource_l, iz, ix, Kokkos::ALL)); + sources.source_array, isource_l, Kokkos::ALL, iz, ix)); // Source time function // For acoustic medium, forward simulation, divide by kappa diff --git a/include/kokkos_abstractions.h b/include/kokkos_abstractions.h index 6cf4368c..69cd2904 100644 --- a/include/kokkos_abstractions.h +++ b/include/kokkos_abstractions.h @@ -572,6 +572,9 @@ template struct array_type { template KOKKOS_INLINE_FUNCTION array_type(const Kokkos::View view) { +#ifndef NDEBUG + assert(view.extent(0) == N); +#endif for (int i = 0; i < N; ++i) { data[i] = view(i); } diff --git a/include/source/adjoint_source.hpp b/include/source/adjoint_source.hpp index db87193e..a95c7c5d 100644 --- a/include/source/adjoint_source.hpp +++ b/include/source/adjoint_source.hpp @@ -28,6 +28,8 @@ class adjoint_source : public source { return specfem::wavefield::type::adjoint; } + std::string print() const override; + private: std::string station_name; std::string network_name; diff --git a/include/source/external.hpp b/include/source/external.hpp new file mode 100644 index 00000000..2e22232f --- /dev/null +++ b/include/source/external.hpp @@ -0,0 +1,39 @@ +#ifndef _SPECFEM_SOURCES_EXTERNAL_HPP1_ +#define _SPECFEM_SOURCES_EXTERNAL_HPP1_ + +#include "compute/compute_mesh.hpp" +#include "compute/compute_partial_derivatives.hpp" +#include "compute/properties/properties.hpp" +#include "source.hpp" +#include "yaml-cpp/yaml.h" + +namespace specfem { +namespace sources { +class external : public source { +public: + external(){}; + + external(YAML::Node &Node, const int nsteps, const type_real dt, + const specfem::wavefield::type wavefield_type) + : wavefield_type(wavefield_type), specfem::sources::source(Node, nsteps, + dt){}; + + void compute_source_array( + const specfem::compute::mesh &mesh, + const specfem::compute::partial_derivatives &partial_derivatives, + const specfem::compute::properties &properties, + specfem::kokkos::HostView3d source_array) override; + + specfem::wavefield::type get_wavefield_type() const override { + return wavefield_type; + } + + std::string print() const override; + +private: + specfem::wavefield::type wavefield_type; +}; +} // namespace sources +} // namespace specfem + +#endif /* _SPECFEM_SOURCES_EXTERNAL_HPP_ */ diff --git a/include/source/interface.hpp b/include/source/interface.hpp index 2de2f2c8..61ac13cd 100644 --- a/include/source/interface.hpp +++ b/include/source/interface.hpp @@ -2,6 +2,7 @@ #define _SOURCE_INTERFACE_HPP #include "adjoint_source.hpp" +#include "external.hpp" #include "force_source.hpp" #include "moment_tensor_source.hpp" #include "read_sources.hpp" diff --git a/include/source_time_function/external.hpp b/include/source_time_function/external.hpp index cefd0b14..11addfb0 100644 --- a/include/source_time_function/external.hpp +++ b/include/source_time_function/external.hpp @@ -34,6 +34,8 @@ class external : public stf { return ss.str(); } + type_real get_t0() const override { return this->__t0; } + private: int __nsteps; type_real __t0; diff --git a/src/reader/seismogram.cpp b/src/reader/seismogram.cpp index fe31caa1..b1cc5881 100644 --- a/src/reader/seismogram.cpp +++ b/src/reader/seismogram.cpp @@ -20,7 +20,7 @@ void specfem::reader::seismogram::read() { while (std::getline(file, line)) { std::istringstream iss(line); - type_real time, value; + double time, value; if (!(iss >> time >> value)) { throw std::runtime_error("Seismogram file " + filename + " is not formatted correctly"); @@ -28,6 +28,7 @@ void specfem::reader::seismogram::read() { source_time_function(nsteps, 0) = time; source_time_function(nsteps, 1) = value; + nsteps++; } file.close(); diff --git a/src/source/adjoint_source.cpp b/src/source/adjoint_source.cpp index 72e686fe..ff43a7d7 100644 --- a/src/source/adjoint_source.cpp +++ b/src/source/adjoint_source.cpp @@ -16,6 +16,7 @@ void specfem::sources::adjoint_source::compute_source_array( const auto N = mesh.quadratures.gll.N; const auto el_type = properties.h_element_types(lcoord.ispec); + const int ncomponents = source_array.extent(0); // Compute lagrange interpolants at the source location auto [hxi_source, hpxi_source] = @@ -31,17 +32,38 @@ void specfem::sources::adjoint_source::compute_source_array( for (int ix = 0; ix < N; ++ix) { hlagrange = hxi_source(ix) * hgamma_source(iz); - if (el_type == specfem::element::medium_tag::acoustic || - (el_type == specfem::element::medium_tag::elastic && - specfem::globals::simulation_wave == specfem::wave::sh)) { + if (el_type == specfem::element::medium_tag::acoustic) { + if (ncomponents != 1) { + throw std::runtime_error( + "Force source requires 1 component for acoustic medium"); + } source_array(0, iz, ix) = hlagrange; - source_array(1, iz, ix) = 0; - } else if ((el_type == specfem::element::medium_tag::elastic && - specfem::globals::simulation_wave == specfem::wave::p_sv) || - el_type == specfem::element::medium_tag::poroelastic) { - source_array(0, iz, ix) = hlagrange; - source_array(1, iz, ix) = hlagrange; + } else if ((el_type == specfem::element::medium_tag::elastic) || + (el_type == specfem::element::medium_tag::poroelastic)) { + if (ncomponents != 2) { + throw std::runtime_error( + "Force source requires 2 components for elastic medium"); + } + if (specfem::globals::simulation_wave == specfem::wave::sh) { + source_array(0, iz, ix) = hlagrange; + source_array(1, iz, ix) = 0; + } else { + source_array(0, iz, ix) = hlagrange; + source_array(1, iz, ix) = hlagrange; + } } } } }; + +std::string specfem::sources::adjoint_source::print() const { + + std::ostringstream message; + message << "- Adjoint Source: \n" + << " Source Location: \n" + << " x = " << type_real(this->x) << "\n" + << " z = " << type_real(this->z) << "\n" + << " Source Time Function: \n" + << this->forcing_function->print() << "\n"; + return message.str(); +} diff --git a/src/source/external.cpp b/src/source/external.cpp new file mode 100644 index 00000000..72f8cde3 --- /dev/null +++ b/src/source/external.cpp @@ -0,0 +1,83 @@ +#include "algorithms/locate_point.hpp" +#include "compute/compute_mesh.hpp" +#include "compute/compute_partial_derivatives.hpp" +#include "compute/properties/properties.hpp" +#include "enumerations/specfem_enums.hpp" +#include "globals.h" +#include "kokkos_abstractions.h" +#include "point/coordinates.hpp" +#include "quadrature/interface.hpp" +#include "source/interface.hpp" +#include "source_time_function/interface.hpp" +#include "specfem_mpi/interface.hpp" +#include "specfem_setup.hpp" +// #include "utilities.cpp" +#include "yaml-cpp/yaml.h" +#include + +void specfem::sources::external::compute_source_array( + const specfem::compute::mesh &mesh, + const specfem::compute::partial_derivatives &partial_derivatives, + const specfem::compute::properties &properties, + specfem::kokkos::HostView3d source_array) { + + specfem::point::gcoord2 coord = specfem::point::gcoord2(this->x, this->z); + auto lcoord = specfem::algorithms::locate_point(coord, mesh); + + const auto xi = mesh.quadratures.gll.h_xi; + const auto gamma = mesh.quadratures.gll.h_xi; + const auto N = mesh.quadratures.gll.N; + + const auto el_type = properties.h_element_types(lcoord.ispec); + const int ncomponents = source_array.extent(0); + + // Compute lagrange interpolants at the source location + auto [hxi_source, hpxi_source] = + specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants( + lcoord.xi, N, xi); + auto [hgamma_source, hpgamma_source] = + specfem::quadrature::gll::Lagrange::compute_lagrange_interpolants( + lcoord.gamma, N, gamma); + + type_real hlagrange; + + for (int iz = 0; iz < N; ++iz) { + for (int ix = 0; ix < N; ++ix) { + hlagrange = hxi_source(ix) * hgamma_source(iz); + + if (el_type == specfem::element::medium_tag::acoustic) { + if (ncomponents != 1) { + throw std::runtime_error( + "Force source requires 1 component for acoustic medium"); + } + source_array(0, iz, ix) = hlagrange; + } else if ((el_type == specfem::element::medium_tag::elastic) || + (el_type == specfem::element::medium_tag::poroelastic)) { + if (ncomponents != 2) { + throw std::runtime_error( + "Force source requires 2 components for elastic medium"); + } + if (specfem::globals::simulation_wave == specfem::wave::sh) { + source_array(0, iz, ix) = hlagrange; + source_array(1, iz, ix) = 0; + } else { + source_array(0, iz, ix) = hlagrange; + source_array(1, iz, ix) = hlagrange; + } + } + } + } +}; + +std::string specfem::sources::external::print() const { + + std::ostringstream message; + message << "- External Source: \n" + << " Source Location: \n" + << " x = " << type_real(this->x) << "\n" + << " z = " << type_real(this->z) << "\n" + << " Source Time Function: \n" + << this->forcing_function->print() << "\n"; + + return message.str(); +} diff --git a/src/source/force_source.cpp b/src/source/force_source.cpp index d5ffffc6..643c4ad2 100644 --- a/src/source/force_source.cpp +++ b/src/source/force_source.cpp @@ -29,7 +29,7 @@ void specfem::sources::force::compute_source_array( const auto N = mesh.quadratures.gll.N; const auto el_type = properties.h_element_types(lcoord.ispec); - const int ncomponents = source_array.extent(2); + const int ncomponents = source_array.extent(0); // Compute lagrange interpolants at the source location auto [hxi_source, hpxi_source] = diff --git a/src/source/moment_tensor_source.cpp b/src/source/moment_tensor_source.cpp index 114b456c..5f4ab33c 100644 --- a/src/source/moment_tensor_source.cpp +++ b/src/source/moment_tensor_source.cpp @@ -33,7 +33,7 @@ void specfem::sources::moment_tensor::compute_source_array( "Moment tensor source not implemented for acoustic medium"); } - const int ncomponents = source_array.extent(2); + const int ncomponents = source_array.extent(0); if (ncomponents != 2) { throw std::runtime_error( "Moment tensor source requires 2 components for elastic medium"); diff --git a/src/source/read_sources.cpp b/src/source/read_sources.cpp index 5017904b..798585fa 100644 --- a/src/source/read_sources.cpp +++ b/src/source/read_sources.cpp @@ -1,3 +1,5 @@ +#include "source/external.hpp" +#include "source/force_source.hpp" #include "source/interface.hpp" #include "source/moment_tensor_source.hpp" #include "specfem_setup.hpp" @@ -39,6 +41,9 @@ specfem::sources::read_sources( } else if (YAML::Node moment_tensor_source = N["moment-tensor"]) { sources.push_back(std::make_shared( moment_tensor_source, nsteps, dt, source_wavefield_type)); + } else if (YAML::Node external = N["user-defined"]) { + sources.push_back(std::make_shared( + external, nsteps, dt, source_wavefield_type)); } else if (YAML::Node adjoint_source = N["adjoint-source"]) { sources.push_back(std::make_shared( adjoint_source, nsteps, dt)); @@ -59,6 +64,7 @@ specfem::sources::read_sources( type_real t0 = std::numeric_limits::max(); for (auto &source : sources) { type_real cur_t0 = source->get_t0(); + std::cout << cur_t0 << std::endl; if (cur_t0 < t0) { t0 = cur_t0; } diff --git a/src/source_time_function/external.cpp b/src/source_time_function/external.cpp index b65b764e..1f30ae9c 100644 --- a/src/source_time_function/external.cpp +++ b/src/source_time_function/external.cpp @@ -10,7 +10,7 @@ specfem::forcing_function::external::external(const YAML::Node &external, const int nsteps, const type_real dt) : __nsteps(nsteps), __dt(dt), - __filename(external["filename"].as()) { + __filename(external["stf-file"].as()) { if ((external["format"].as() == "ascii") || (external["format"].as() == "ASCII")) { @@ -19,7 +19,9 @@ specfem::forcing_function::external::external(const YAML::Node &external, throw std::runtime_error("Only ASCII format is supported"); } - std::vector extentions = { "BXX.semd", "BXY.semd" }; + std::vector extentions = { ".BXX.semd", ".BXY.semd" }; + + int file_components = 0; // get t0 from file for (int icomp = 0; icomp < 2; ++icomp) { @@ -35,10 +37,17 @@ specfem::forcing_function::external::external(const YAML::Node &external, " is not formatted correctly"); } this->__t0 = time; + file_components++; } file.close(); } + if (file_components != 1) { + throw std::runtime_error( + "Wrong number for files found for external source: " + + this->__filename); + } + return; } @@ -50,9 +59,9 @@ void specfem::forcing_function::external::compute_source_time_function( const auto extention = [&ncomponents]() -> std::vector { if (ncomponents == 2) { - return { "BXX.semd", "BXZ.semd" }; + return { ".BXX.semd", ".BXZ.semd" }; } else if (ncomponents == 1) { - return { "BXY.semd" }; + return { ".BXY.semd" }; } else { throw std::runtime_error("Invalid number of components"); } @@ -69,7 +78,7 @@ void specfem::forcing_function::external::compute_source_time_function( } for (int icomp = 0; icomp < ncomponents; ++icomp) { - specfem::kokkos::HostView2d data("external", nsteps); + specfem::kokkos::HostView2d data("external", nsteps, 2); specfem::reader::seismogram reader( this->__filename + extention[icomp], specfem::enums::seismogram::format::ascii, data);