From fe6635e0e3cf13632b4f8ffdbd133d5440d8e331 Mon Sep 17 00:00:00 2001 From: Rohit Kakodkar Date: Wed, 19 Feb 2025 15:59:59 -0500 Subject: [PATCH] Updates mesh reader for shear waves - Implements point::properties for shear waves. - To do this I create a base class with medium_tag = specfem::element::medium_tag::elastic - Implements medium::properties for shear waves - Again this is implemented with a base class with medium_tag = specfem::element::medium_tag::elastic - Update specfem::IO::read_mesh signature to include a elastic_wave argument. The elastic_wave would be set up using specfem_config.yaml - Updates tests for the mesher --- include/IO/interface.hpp | 3 +- .../impl/fortran/read_material_properties.hpp | 1 + include/enumerations/medium.hpp | 8 +- include/enumerations/specfem_enums.hpp | 2 + .../dim2/elastic/anisotropic/properties.hpp | 60 +++++- .../anisotropic/properties_container.hpp | 8 +- .../dim2/elastic/isotropic/properties.hpp | 54 +++++- .../isotropic/properties_container.hpp | 8 +- include/mesh/materials/materials.hpp | 13 ++ include/point/properties.hpp | 174 +++++++++++++++++- src/IO/mesh.cpp | 5 +- .../impl/fortran/read_material_properties.cpp | 110 ++++++++--- src/execute.cpp | 3 +- src/mesh/materials/materials.cpp | 19 ++ .../algorithms/interpolate_function.cpp | 3 +- tests/unit-tests/algorithms/locate.cpp | 3 +- .../assembly/properties/properties.cpp | 70 ++++++- .../assembly/test_fixture/test_fixture.cpp | 3 +- .../compute/index/compute_tests.cpp | 4 +- .../compute_partial_derivatives_tests.cpp | 4 +- .../Newmark/newmark_tests.cpp | 3 +- .../unit-tests/domain/rmass_inverse_tests.cpp | 3 +- tests/unit-tests/mesh/materials/materials.cpp | 54 +++++- tests/unit-tests/mesh/test_config.yaml | 49 +++-- .../mesh/test_fixture/test_fixture.cpp | 4 +- .../mesh/test_fixture/test_fixture.hpp | 37 ++++ 26 files changed, 618 insertions(+), 87 deletions(-) diff --git a/include/IO/interface.hpp b/include/IO/interface.hpp index 8c1429fb5..6441e1792 100644 --- a/include/IO/interface.hpp +++ b/include/IO/interface.hpp @@ -22,7 +22,8 @@ namespace IO { * */ specfem::mesh::mesh -read_mesh(const std::string filename, const specfem::MPI::MPI *mpi); +read_mesh(const std::string filename, const specfem::enums::elastic_wave wave, + const specfem::MPI::MPI *mpi); /** * @brief Read station file diff --git a/include/IO/mesh/impl/fortran/read_material_properties.hpp b/include/IO/mesh/impl/fortran/read_material_properties.hpp index ecf64b068..2efd3305f 100644 --- a/include/IO/mesh/impl/fortran/read_material_properties.hpp +++ b/include/IO/mesh/impl/fortran/read_material_properties.hpp @@ -26,6 +26,7 @@ namespace fortran { specfem::mesh::materials read_material_properties( std::ifstream &stream, const int numat, const int nspec, + const specfem::enums::elastic_wave wave, const specfem::kokkos::HostView2d knods, const specfem::MPI::MPI *mpi); } // namespace fortran diff --git a/include/enumerations/medium.hpp b/include/enumerations/medium.hpp index 6ef64829a..85ad76cd6 100644 --- a/include/enumerations/medium.hpp +++ b/include/enumerations/medium.hpp @@ -13,7 +13,13 @@ constexpr int ntypes = 3; ///< Number of element types * @brief Medium tag enumeration * */ -enum class medium_tag { elastic_sv, elastic_sh, acoustic, poroelastic }; +enum class medium_tag { + elastic_sv, + elastic_sh, + acoustic, + elastic, + poroelastic +}; /** * @brief Property tag enumeration diff --git a/include/enumerations/specfem_enums.hpp b/include/enumerations/specfem_enums.hpp index 497f748cf..be1ba1877 100644 --- a/include/enumerations/specfem_enums.hpp +++ b/include/enumerations/specfem_enums.hpp @@ -10,6 +10,8 @@ namespace specfem { */ namespace enums { +enum class elastic_wave { p_sv, sh }; + /** * @brief Cartesian axes enumeration * diff --git a/include/medium/dim2/elastic/anisotropic/properties.hpp b/include/medium/dim2/elastic/anisotropic/properties.hpp index 2cdb1d43b..587ab4a7d 100644 --- a/include/medium/dim2/elastic/anisotropic/properties.hpp +++ b/include/medium/dim2/elastic/anisotropic/properties.hpp @@ -17,7 +17,7 @@ namespace medium { * */ template <> -class properties { public: constexpr static auto dimension = @@ -77,7 +77,7 @@ class properties &other) const { return (std::abs(this->density - other.density) < 1e-6 && @@ -100,7 +100,7 @@ class properties &other) const { return !(*this == other); @@ -151,5 +151,59 @@ class properties +class properties + : public properties { + +public: + constexpr static auto dimension = + specfem::dimension::type::dim2; ///< Dimension of the material + constexpr static auto medium_tag = + specfem::element::medium_tag::elastic_sv; ///< Medium tag + constexpr static auto property_tag = + specfem::element::property_tag::anisotropic; ///< Property tag + + properties(const type_real &density, const type_real &c11, + const type_real &c13, const type_real &c15, const type_real &c33, + const type_real &c35, const type_real &c55, const type_real &c12, + const type_real &c23, const type_real &c25, + const type_real &Qkappa, const type_real &Qmu) + : properties( + density, c11, c13, c15, c33, c35, c55, c12, c23, c25, Qkappa, Qmu) { + } + + properties() = default; +}; + +template <> +class properties + : public properties { + +public: + constexpr static auto dimension = + specfem::dimension::type::dim2; ///< Dimension of the material + constexpr static auto medium_tag = + specfem::element::medium_tag::elastic_sh; ///< Medium tag + constexpr static auto property_tag = + specfem::element::property_tag::anisotropic; ///< Property tag + + properties(const type_real &density, const type_real &c11, + const type_real &c13, const type_real &c15, const type_real &c33, + const type_real &c35, const type_real &c55, const type_real &c12, + const type_real &c23, const type_real &c25, + const type_real &Qkappa, const type_real &Qmu) + : properties( + density, c11, c13, c15, c33, c35, c55, c12, c23, c25, Qkappa, Qmu) { + } + + properties() = default; +}; + } // namespace medium } // namespace specfem diff --git a/include/medium/dim2/elastic/anisotropic/properties_container.hpp b/include/medium/dim2/elastic/anisotropic/properties_container.hpp index 97f5507a2..3f3531ffc 100644 --- a/include/medium/dim2/elastic/anisotropic/properties_container.hpp +++ b/include/medium/dim2/elastic/anisotropic/properties_container.hpp @@ -265,8 +265,8 @@ struct properties_container -class properties { public: constexpr static auto dimension = specfem::dimension::type::dim2; ///< Dimension of the material constexpr static auto medium_tag = - specfem::element::medium_tag::elastic_sv; ///< Medium tag + specfem::element::medium_tag::elastic; ///< Medium tag constexpr static auto property_tag = specfem::element::property_tag::isotropic; ///< Property tag @@ -73,7 +73,7 @@ class properties &other) const { @@ -91,7 +91,7 @@ class properties &other) const { return !(*this == other); @@ -142,5 +142,51 @@ class properties +class properties + : public properties { +public: + constexpr static auto dimension = + specfem::dimension::type::dim2; ///< Dimension of the material + constexpr static auto medium_tag = + specfem::element::medium_tag::elastic_sv; ///< Medium tag + constexpr static auto property_tag = + specfem::element::property_tag::isotropic; ///< Property tag + + properties(const type_real &density, const type_real &cs, const type_real &cp, + const type_real &Qkappa, const type_real &Qmu, + const type_real &compaction_grad) + : properties( + density, cs, cp, Qkappa, Qmu, compaction_grad) {} + + properties() = default; +}; + +template <> +class properties + : public properties { +public: + constexpr static auto dimension = + specfem::dimension::type::dim2; ///< Dimension of the material + constexpr static auto medium_tag = + specfem::element::medium_tag::elastic_sh; ///< Medium tag + constexpr static auto property_tag = + specfem::element::property_tag::isotropic; ///< Property tag + + properties(const type_real &density, const type_real &cs, const type_real &cp, + const type_real &Qkappa, const type_real &Qmu, + const type_real &compaction_grad) + : properties( + density, cs, cp, Qkappa, Qmu, compaction_grad) {} + + properties() = default; +}; + } // namespace medium } // namespace specfem diff --git a/include/medium/dim2/elastic/isotropic/properties_container.hpp b/include/medium/dim2/elastic/isotropic/properties_container.hpp index cf1dc2d96..4537f5ab9 100644 --- a/include/medium/dim2/elastic/isotropic/properties_container.hpp +++ b/include/medium/dim2/elastic/isotropic/properties_container.hpp @@ -185,8 +185,8 @@ struct properties_container elastic_sv_anisotropic; ///< Elastic anisotropic material properties + specfem::mesh::materials::material + elastic_sh_isotropic; ///< Elastic isotropic material properties + + specfem::mesh::materials::material< + specfem::element::medium_tag::elastic_sh, + specfem::element::property_tag::anisotropic> + elastic_sh_anisotropic; ///< Elastic anisotropic material properties + specfem::mesh::materials::material acoustic_isotropic; ///< Acoustic isotropic material properties @@ -104,6 +113,10 @@ struct materials { specfem::element::property_tag::isotropic>, specfem::medium::material, + specfem::medium::material, + specfem::medium::material, specfem::medium::material > operator[](const int index) const; diff --git a/include/point/properties.hpp b/include/point/properties.hpp index 0c3efa65f..08fb6b39c 100644 --- a/include/point/properties.hpp +++ b/include/point/properties.hpp @@ -27,7 +27,7 @@ struct properties; */ template struct properties { /** @@ -132,7 +132,7 @@ struct properties struct properties { /** @@ -200,7 +200,7 @@ struct properties{}) {} @@ -240,6 +240,174 @@ struct properties +struct properties + : public properties { + + /** + * @name Typedefs + * + */ + ///@{ + constexpr static bool is_point_properties = true; + using simd = specfem::datatype::simd; ///< SIMD type + constexpr static auto dimension = specfem::dimension::type::dim2; + constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv; + constexpr static auto property_tag = + specfem::element::property_tag::isotropic; + using value_type = + typename simd::datatype; ///< Value type to store properties + ///@} + + KOKKOS_FUNCTION + properties() = default; + + KOKKOS_FUNCTION + properties(const value_type &lambdaplus2mu, const value_type &mu, + const value_type &rho) + : properties( + lambdaplus2mu, mu, rho) {} + + properties(const value_type value) + : properties(value) {} +}; + +template +struct properties + : public properties { + + /** + * @name Typedefs + * + */ + ///@{ + constexpr static bool is_point_properties = true; + using simd = specfem::datatype::simd; ///< SIMD type + constexpr static auto dimension = specfem::dimension::type::dim2; + constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sv; + constexpr static auto property_tag = + specfem::element::property_tag::anisotropic; + using value_type = + typename simd::datatype; ///< Value type to store properties + ///@} + + KOKKOS_FUNCTION + properties() = default; + + KOKKOS_FUNCTION + properties(const value_type &c11, const value_type &c13, + const value_type &c15, const value_type &c33, + const value_type &c35, const value_type &c55, + const value_type &c12, const value_type &c23, + const value_type &c25, const value_type &rho) + : properties( + c11, c13, c15, c33, c35, c55, c12, c23, c25, rho) {} + + properties(const value_type value) + : properties( + value) {} +}; + +template +struct properties + : public properties { + + /** + * @name Typedefs + * + */ + ///@{ + constexpr static bool is_point_properties = true; + using simd = specfem::datatype::simd; ///< SIMD type + constexpr static auto dimension = specfem::dimension::type::dim2; + constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh; + constexpr static auto property_tag = + specfem::element::property_tag::isotropic; + using value_type = + typename simd::datatype; ///< Value type to store properties + ///@} + + KOKKOS_FUNCTION + properties() = default; + + KOKKOS_FUNCTION + properties(const value_type &lambdaplus2mu, const value_type &mu, + const value_type &rho) + : properties( + lambdaplus2mu, mu, rho) {} + + properties(const value_type value) + : properties(value) {} +}; + +template +struct properties + : public properties { + + /** + * @name Typedefs + * + */ + ///@{ + constexpr static bool is_point_properties = true; + using simd = specfem::datatype::simd; ///< SIMD type + constexpr static auto dimension = specfem::dimension::type::dim2; + constexpr static auto medium_tag = specfem::element::medium_tag::elastic_sh; + constexpr static auto property_tag = + specfem::element::property_tag::anisotropic; + using value_type = + typename simd::datatype; ///< Value type to store properties + ///@} + + KOKKOS_FUNCTION + properties() = default; + + KOKKOS_FUNCTION + properties(const value_type &c11, const value_type &c13, + const value_type &c15, const value_type &c33, + const value_type &c35, const value_type &c55, + const value_type &c12, const value_type &c23, + const value_type &c25, const value_type &rho) + : properties( + c11, c13, c15, c33, c35, c55, c12, c23, c25, rho) {} + + properties(const value_type value) + : properties( + value) {} +}; + /** * @brief Template specialization for 2D isotropic acoustic media * diff --git a/src/IO/mesh.cpp b/src/IO/mesh.cpp index d57a143c3..3ad22eb54 100644 --- a/src/IO/mesh.cpp +++ b/src/IO/mesh.cpp @@ -25,6 +25,7 @@ specfem::mesh::mesh specfem::IO::read_mesh(const std::string filename, + const specfem::enums::elastic_wave wave, const specfem::MPI::MPI *mpi) { // Declaring empty mesh objects @@ -87,8 +88,8 @@ specfem::IO::read_mesh(const std::string filename, try { mesh.materials = specfem::IO::mesh::impl::fortran::read_material_properties( - stream, mesh.parameters.numat, mesh.nspec, mesh.control_nodes.knods, - mpi); + stream, mesh.parameters.numat, mesh.nspec, wave, + mesh.control_nodes.knods, mpi); } catch (std::runtime_error &e) { throw; } diff --git a/src/IO/mesh/impl/fortran/read_material_properties.cpp b/src/IO/mesh/impl/fortran/read_material_properties.cpp index 8e43481b6..843c75e7e 100644 --- a/src/IO/mesh/impl/fortran/read_material_properties.cpp +++ b/src/IO/mesh/impl/fortran/read_material_properties.cpp @@ -10,6 +10,7 @@ // Define some constants for the material properties constexpr auto elastic_sv = specfem::element::medium_tag::elastic_sv; +constexpr auto elastic_sh = specfem::element::medium_tag::elastic_sh; constexpr auto isotropic = specfem::element::property_tag::isotropic; constexpr auto anisotropic = specfem::element::property_tag::anisotropic; constexpr auto acoustic = specfem::element::medium_tag::acoustic; @@ -23,13 +24,27 @@ struct input_holder { std::vector read_materials( std::ifstream &stream, const int numat, + const specfem::enums::elastic_wave wave, specfem::mesh::materials::material - &elastic_isotropic, + &elastic_sv_isotropic, + specfem::mesh::materials::material + &elastic_sh_isotropic, specfem::mesh::materials::material &acoustic_isotropic, specfem::mesh::materials::material - &elastic_anisotropic, + &elastic_sv_anisotropic, + specfem::mesh::materials::material + &elastic_sh_anisotropic, const specfem::MPI::MPI *mpi) { + const specfem::element::medium_tag elastic = [wave]() { + if (wave == specfem::enums::elastic_wave::p_sv) + return specfem::element::medium_tag::elastic_sv; + else if (wave == specfem::enums::elastic_wave::sh) + return specfem::element::medium_tag::elastic_sh; + else + throw std::runtime_error("Elastic wave type not supported"); + }(); + input_holder read_values; std::ostringstream message; @@ -47,17 +62,25 @@ std::vector read_materials( // Section for elastic isotropic std::vector > - l_elastic_isotropic; + l_elastic_sv_isotropic; + + std::vector > + l_elastic_sh_isotropic; - l_elastic_isotropic.reserve(numat); + l_elastic_sv_isotropic.reserve(numat); + l_elastic_sh_isotropic.reserve(numat); int index_elastic_isotropic = 0; // Section for elastic anisotropic std::vector > - l_elastic_anisotropic; + l_elastic_sv_anisotropic; + + std::vector > + l_elastic_sh_anisotropic; - l_elastic_anisotropic.reserve(numat); + l_elastic_sv_anisotropic.reserve(numat); + l_elastic_sh_anisotropic.reserve(numat); int index_elastic_anisotropic = 0; @@ -121,18 +144,25 @@ std::vector read_materials( const type_real Qkappa = read_values.val5; const type_real Qmu = read_values.val6; - specfem::medium::material - elastic_isotropic_holder(density, cs, cp, Qkappa, Qmu, - compaction_grad); + if (wave == specfem::enums::elastic_wave::p_sv) { + specfem::medium::material + elastic_isotropic_holder(density, cs, cp, Qkappa, Qmu, + compaction_grad); - elastic_isotropic_holder.print(); + elastic_isotropic_holder.print(); + l_elastic_sv_isotropic.push_back(elastic_isotropic_holder); + } else { + specfem::medium::material + elastic_isotropic_holder(density, cs, cp, Qkappa, Qmu, + compaction_grad); - l_elastic_isotropic.push_back(elastic_isotropic_holder); + elastic_isotropic_holder.print(); + l_elastic_sh_isotropic.push_back(elastic_isotropic_holder); + } index_mapping[i] = specfem::mesh::materials::material_specification( - specfem::element::medium_tag::elastic_sv, - specfem::element::property_tag::isotropic, index_elastic_isotropic, - read_values.n - 1); + elastic, specfem::element::property_tag::isotropic, + index_elastic_isotropic, read_values.n - 1); index_elastic_isotropic++; } @@ -152,17 +182,26 @@ std::vector read_materials( const type_real Qkappa = read_values.val11; const type_real Qmu = read_values.val12; - specfem::medium::material - elastic_anisotropic_holder(density, c11, c13, c15, c33, c35, c55, c12, - c23, c25, Qkappa, Qmu); + if (wave == specfem::enums::elastic_wave::p_sv) { + + specfem::medium::material + elastic_anisotropic_holder(density, c11, c13, c15, c33, c35, c55, + c12, c23, c25, Qkappa, Qmu); - elastic_anisotropic_holder.print(); + elastic_anisotropic_holder.print(); + l_elastic_sv_anisotropic.push_back(elastic_anisotropic_holder); + } else { + + specfem::medium::material + elastic_anisotropic_holder(density, c11, c13, c15, c33, c35, c55, + c12, c23, c25, Qkappa, Qmu); - l_elastic_anisotropic.push_back(elastic_anisotropic_holder); + elastic_anisotropic_holder.print(); + l_elastic_sh_anisotropic.push_back(elastic_anisotropic_holder); + } index_mapping[i] = specfem::mesh::materials::material_specification( - specfem::element::medium_tag::elastic_sv, - specfem::element::property_tag::anisotropic, + elastic, specfem::element::property_tag::anisotropic, index_elastic_anisotropic, read_values.n - 1); index_elastic_anisotropic++; @@ -172,16 +211,25 @@ std::vector read_materials( } } - assert(l_elastic_isotropic.size() + l_acoustic_isotropic.size() + - l_elastic_anisotropic.size() == - numat); + assert((l_elastic_sv_isotropic.size() + l_elastic_sh_isotropic.size() + + l_acoustic_isotropic.size() + l_elastic_sv_anisotropic.size() + + l_elastic_sh_anisotropic.size()) == numat); - elastic_isotropic = specfem::mesh::materials::material( - l_elastic_isotropic.size(), l_elastic_isotropic); + elastic_sv_isotropic = + specfem::mesh::materials::material( + l_elastic_sv_isotropic.size(), l_elastic_sv_isotropic); - elastic_anisotropic = + elastic_sh_isotropic = + specfem::mesh::materials::material( + l_elastic_sh_isotropic.size(), l_elastic_sh_isotropic); + + elastic_sv_anisotropic = specfem::mesh::materials::material( - l_elastic_anisotropic.size(), l_elastic_anisotropic); + l_elastic_sv_anisotropic.size(), l_elastic_sv_anisotropic); + + elastic_sh_anisotropic = + specfem::mesh::materials::material( + l_elastic_sh_anisotropic.size(), l_elastic_sh_anisotropic); acoustic_isotropic = specfem::mesh::materials::material( l_acoustic_isotropic.size(), l_acoustic_isotropic); @@ -234,6 +282,7 @@ void read_material_indices( specfem::mesh::materials specfem::IO::mesh::impl::fortran::read_material_properties( std::ifstream &stream, const int numat, const int nspec, + const specfem::enums::elastic_wave wave, const specfem::kokkos::HostView2d knods, const specfem::MPI::MPI *mpi) { @@ -242,8 +291,9 @@ specfem::IO::mesh::impl::fortran::read_material_properties( // Read material properties auto index_mapping = read_materials( - stream, numat, materials.elastic_sv_isotropic, - materials.acoustic_isotropic, materials.elastic_sv_anisotropic, mpi); + stream, numat, wave, materials.elastic_sv_isotropic, + materials.elastic_sh_isotropic, materials.acoustic_isotropic, + materials.elastic_sv_anisotropic, materials.elastic_sh_anisotropic, mpi); // Read material indices read_material_indices(stream, nspec, numat, index_mapping, diff --git a/src/execute.cpp b/src/execute.cpp index 138d51902..6418fb2fc 100644 --- a/src/execute.cpp +++ b/src/execute.cpp @@ -42,7 +42,8 @@ void execute( // Read mesh and materials // -------------------------------------------------------------- const auto quadrature = setup.instantiate_quadrature(); - const auto mesh = specfem::IO::read_mesh(database_filename, mpi); + const auto mesh = specfem::IO::read_mesh( + database_filename, specfem::enums::elastic_wave::p_sv, mpi); // -------------------------------------------------------------- // -------------------------------------------------------------- diff --git a/src/mesh/materials/materials.cpp b/src/mesh/materials/materials.cpp index 346e40aa7..6aedbcdef 100644 --- a/src/mesh/materials/materials.cpp +++ b/src/mesh/materials/materials.cpp @@ -9,6 +9,10 @@ std::variant< specfem::element::property_tag::isotropic>, specfem::medium::material, + specfem::medium::material, + specfem::medium::material, specfem::medium::material > specfem::mesh::materials::operator[](const int index) const { @@ -32,6 +36,21 @@ specfem::mesh::materials::operator[](const int index) const { return this->elastic_sv_anisotropic .material_properties[material_specification.index]; + } else if (material_specification.type == + specfem::element::medium_tag::elastic_sh && + material_specification.property == + specfem::element::property_tag::isotropic) { + return this->elastic_sh_isotropic + .material_properties[material_specification.index]; + + // Return elastic_sh anisotropic + } else if (material_specification.type == + specfem::element::medium_tag::elastic_sh && + material_specification.property == + specfem::element::property_tag::anisotropic) { + return this->elastic_sh_anisotropic + .material_properties[material_specification.index]; + // Return acoustic isotropic } else if (material_specification.type == specfem::element::medium_tag::acoustic && diff --git a/tests/unit-tests/algorithms/interpolate_function.cpp b/tests/unit-tests/algorithms/interpolate_function.cpp index 3bfecefa6..3dbac593c 100644 --- a/tests/unit-tests/algorithms/interpolate_function.cpp +++ b/tests/unit-tests/algorithms/interpolate_function.cpp @@ -21,7 +21,8 @@ TEST(ALGORITHMS, interpolate_function) { // Read Mesh database specfem::MPI::MPI *mpi = MPIEnvironment::get_mpi(); - specfem::mesh::mesh mesh = specfem::IO::read_mesh(database_file, mpi); + specfem::mesh::mesh mesh = specfem::IO::read_mesh( + database_file, specfem::enums::elastic_wave::p_sv, mpi); constexpr int N = 5; diff --git a/tests/unit-tests/algorithms/locate.cpp b/tests/unit-tests/algorithms/locate.cpp index 607732c16..83e91faf0 100644 --- a/tests/unit-tests/algorithms/locate.cpp +++ b/tests/unit-tests/algorithms/locate.cpp @@ -14,7 +14,8 @@ TEST(ALGORITHMS, locate_point) { // Read Mesh database specfem::MPI::MPI *mpi = MPIEnvironment::get_mpi(); - specfem::mesh::mesh mesh = specfem::IO::read_mesh(database_file, mpi); + specfem::mesh::mesh mesh = specfem::IO::read_mesh( + database_file, specfem::enums::elastic_wave::p_sv, mpi); // Quadratures specfem::quadrature::gll::gll gll(0.0, 0.0, 5); diff --git a/tests/unit-tests/assembly/properties/properties.cpp b/tests/unit-tests/assembly/properties/properties.cpp index 7be7c076c..242c309ef 100644 --- a/tests/unit-tests/assembly/properties/properties.cpp +++ b/tests/unit-tests/assembly/properties/properties.cpp @@ -31,10 +31,9 @@ std::string get_error_message( template <> std::string get_error_message( - const specfem::point::properties &point_property, + const specfem::point::properties< + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, + specfem::element::property_tag::isotropic, false> &point_property, const type_real value, const int mode) { std::ostringstream message; @@ -49,8 +48,7 @@ std::string get_error_message( template <> std::string get_error_message( const specfem::point::properties< - specfem::dimension::type::dim2, - specfem::element::medium_tag::elastic_sv, + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, specfem::element::property_tag::anisotropic, false> &point_property, const type_real value, const int mode) { std::ostringstream message; @@ -67,6 +65,66 @@ std::string get_error_message( return message.str(); } +template <> +std::string get_error_message( + const specfem::point::properties &point_property, + const type_real value, const int mode) { + + return get_error_message( + static_cast >(point_property), + value, mode); +} + +template <> +std::string get_error_message( + const specfem::point::properties< + specfem::dimension::type::dim2, + specfem::element::medium_tag::elastic_sv, + specfem::element::property_tag::anisotropic, false> &point_property, + const type_real value, const int mode) { + + return get_error_message( + static_cast >(point_property), + value, mode); +} + +template <> +std::string get_error_message( + const specfem::point::properties &point_property, + const type_real value, const int mode) { + + return get_error_message( + static_cast >(point_property), + value, mode); +} + +template <> +std::string get_error_message( + const specfem::point::properties< + specfem::dimension::type::dim2, + specfem::element::medium_tag::elastic_sh, + specfem::element::property_tag::anisotropic, false> &point_property, + const type_real value, const int mode) { + + return get_error_message( + static_cast >(point_property), + value, mode); +} + template <> std::string get_error_message( const specfem::point::properties< diff --git a/tests/unit-tests/assembly/test_fixture/test_fixture.cpp b/tests/unit-tests/assembly/test_fixture/test_fixture.cpp index a338c7f8a..4ff4889c6 100644 --- a/tests/unit-tests/assembly/test_fixture/test_fixture.cpp +++ b/tests/unit-tests/assembly/test_fixture/test_fixture.cpp @@ -32,7 +32,8 @@ ASSEMBLY::ASSEMBLY() { for (auto &Test : Tests) { const auto [database_file, sources_file, stations_file] = Test.get_databases(); - const auto mesh = specfem::IO::read_mesh(database_file, mpi); + const auto mesh = specfem::IO::read_mesh( + database_file, specfem::enums::elastic_wave::p_sv, mpi); this->Meshes.push_back(mesh); diff --git a/tests/unit-tests/compute/index/compute_tests.cpp b/tests/unit-tests/compute/index/compute_tests.cpp index a624d0e10..833ba5d5f 100644 --- a/tests/unit-tests/compute/index/compute_tests.cpp +++ b/tests/unit-tests/compute/index/compute_tests.cpp @@ -81,8 +81,8 @@ TEST(COMPUTE_TESTS, compute_ibool) { specfem::quadrature::quadratures quadratures(gll); // Read mesh generated MESHFEM - specfem::mesh::mesh mesh = - specfem::IO::read_mesh(test_config.database_filename, mpi); + specfem::mesh::mesh mesh = specfem::IO::read_mesh( + test_config.database_filename, specfem::enums::elastic_wave::p_sv, mpi); // Setup compute structs specfem::compute::mesh assembly(mesh.tags, mesh.control_nodes, diff --git a/tests/unit-tests/compute/partial_derivatives/compute_partial_derivatives_tests.cpp b/tests/unit-tests/compute/partial_derivatives/compute_partial_derivatives_tests.cpp index c7bd1a150..3f558e73f 100644 --- a/tests/unit-tests/compute/partial_derivatives/compute_partial_derivatives_tests.cpp +++ b/tests/unit-tests/compute/partial_derivatives/compute_partial_derivatives_tests.cpp @@ -74,8 +74,8 @@ TEST(COMPUTE_TESTS, compute_partial_derivatives) { specfem::quadrature::gll::gll gll(0.0, 0.0, 5); specfem::quadrature::quadratures quadratures(gll); - specfem::mesh::mesh mesh = - specfem::IO::read_mesh(test_config.database_filename, mpi); + specfem::mesh::mesh mesh = specfem::IO::read_mesh( + test_config.database_filename, specfem::enums::elastic_wave::p_sv, mpi); specfem::compute::mesh compute_mesh(mesh.tags, mesh.control_nodes, quadratures); diff --git a/tests/unit-tests/displacement_tests/Newmark/newmark_tests.cpp b/tests/unit-tests/displacement_tests/Newmark/newmark_tests.cpp index d1bfcede0..009b035f3 100644 --- a/tests/unit-tests/displacement_tests/Newmark/newmark_tests.cpp +++ b/tests/unit-tests/displacement_tests/Newmark/newmark_tests.cpp @@ -175,7 +175,8 @@ TEST(DISPLACEMENT_TESTS, newmark_scheme_tests) { const auto quadratures = setup.instantiate_quadrature(); // Read mesh generated MESHFEM - specfem::mesh::mesh mesh = specfem::IO::read_mesh(database_file, mpi); + specfem::mesh::mesh mesh = specfem::IO::read_mesh( + database_file, specfem::enums::elastic_wave::p_sv, mpi); const type_real dt = setup.get_dt(); const int nsteps = setup.get_nsteps(); diff --git a/tests/unit-tests/domain/rmass_inverse_tests.cpp b/tests/unit-tests/domain/rmass_inverse_tests.cpp index 6a4128ab6..880c4b71c 100644 --- a/tests/unit-tests/domain/rmass_inverse_tests.cpp +++ b/tests/unit-tests/domain/rmass_inverse_tests.cpp @@ -117,7 +117,8 @@ TEST(DOMAIN_TESTS, rmass_inverse) { std::cout << "Reading mesh file: " << database_file << std::endl; // Read mesh generated MESHFEM - specfem::mesh::mesh mesh = specfem::IO::read_mesh(database_file, mpi); + specfem::mesh::mesh mesh = specfem::IO::read_mesh( + database_file, specfem::enums::elastic_wave::p_sv, mpi); std::cout << "Setting up sources and receivers" << std::endl; diff --git a/tests/unit-tests/mesh/materials/materials.cpp b/tests/unit-tests/mesh/materials/materials.cpp index 7f6d468b4..331dd53f0 100644 --- a/tests/unit-tests/mesh/materials/materials.cpp +++ b/tests/unit-tests/mesh/materials/materials.cpp @@ -10,22 +10,31 @@ using MaterialVectorType = std::vector, specfem::medium::material, + specfem::medium::material, specfem::medium::material, + specfem::medium::material > >; const static std::unordered_map ground_truth = { - { "Test 1: Simple mesh with flat topography", + { "Test 1: Simple mesh with flat topography (P_SV wave)", MaterialVectorType({ specfem::medium::material< specfem::element::medium_tag::elastic_sv, specfem::element::property_tag::isotropic>(2700.0, 1732.051, 3000.0, 9999, 9999, 0.0) }) }, - { "Test 2: Simple mesh with curved topography", + { "Test 2: Simple mesh with flat topography (SH wave)", + MaterialVectorType({ specfem::medium::material< + specfem::element::medium_tag::elastic_sh, + specfem::element::property_tag::isotropic>(2700.0, 1732.051, 3000.0, + 9999, 9999, 0.0) }) }, + { "Test 3: Simple mesh with curved topography", MaterialVectorType({ specfem::medium::material< specfem::element::medium_tag::elastic_sv, specfem::element::property_tag::isotropic>(2700.0, 1732.051, 3000.0, 9999, 9999, 0.0) }) }, - { "Test 3: Simple mesh with flat ocean bottom", + { "Test 4: Simple mesh with flat ocean bottom", MaterialVectorType({ specfem::medium::material< specfem::element::medium_tag::elastic_sv, specfem::element::property_tag::isotropic>( @@ -36,7 +45,7 @@ const static std::unordered_map 1020.0, 1500, 9999, 9999, 0.0) }) }, - { "Test 4: Simple mesh with curved ocean bottom", + { "Test 5: Simple mesh with curved ocean bottom", MaterialVectorType({ specfem::medium::material< specfem::element::medium_tag::elastic_sv, specfem::element::property_tag::isotropic>( @@ -47,7 +56,7 @@ const static std::unordered_map 1020.0, 1500, 9999, 9999, 0.0) }) }, - { "Test 5: Gmesh Example", + { "Test 6: Gmesh Example", MaterialVectorType({ specfem::medium::material< specfem::element::medium_tag::acoustic, specfem::element::property_tag::isotropic>( @@ -56,11 +65,18 @@ const static std::unordered_map specfem::element::medium_tag::acoustic, specfem::element::property_tag::isotropic>( 1000.0, 1477.0, 10.0, 10.0, 0.0) }) }, - { "Test 6: Homogeneous Elastic Anisotropic Material", + { "Test 7: Homogeneous Elastic Anisotropic Material (P_SV wave)", MaterialVectorType({ specfem::medium::material< specfem::element::medium_tag::elastic_sv, specfem::element::property_tag::anisotropic>( 2700.0, 24299994600.5, 8099996400.35, 0.0, 24299994600.5, 0.0, + 8100001799.8227, 8099996400.35, 8099996400.35, 0.0, 9999, + 9999) }) }, + { "Test 8: Homogeneous Elastic Anisotropic Material (SH wave)", + MaterialVectorType({ specfem::medium::material< + specfem::element::medium_tag::elastic_sh, + specfem::element::property_tag::anisotropic>( + 2700.0, 24299994600.5, 8099996400.35, 0.0, 24299994600.5, 0.0, 8100001799.8227, 8099996400.35, 8099996400.35, 0.0, 9999, 9999) }) } }; @@ -112,6 +128,19 @@ void check_test(const specfem::mesh::materials &computed, error_message << "Material " << index << " is not the same"; throw std::runtime_error(error_message.str()); } + } else if ((type == specfem::element::medium_tag::elastic_sh) && + (property == specfem::element::property_tag::isotropic)) { + const auto icomputed = std::get >(computed[ispec]); + const auto iexpected = std::get >(expected[imaterial]); + if (icomputed != iexpected) { + std::ostringstream error_message; + error_message << "Material " << index << " is not the same"; + throw std::runtime_error(error_message.str()); + } } else if ((type == specfem::element::medium_tag::elastic_sv) && (property == specfem::element::property_tag::anisotropic)) { const auto icomputed = std::get >(computed[ispec]); + const auto iexpected = std::get >(expected[imaterial]); + if (icomputed != iexpected) { + std::ostringstream error_message; + error_message << "Material " << index << " is not the same"; + throw std::runtime_error(error_message.str()); + } } else { throw std::runtime_error("Material type not supported"); } diff --git a/tests/unit-tests/mesh/test_config.yaml b/tests/unit-tests/mesh/test_config.yaml index 07b7b8e5c..3c4853e45 100644 --- a/tests/unit-tests/mesh/test_config.yaml +++ b/tests/unit-tests/mesh/test_config.yaml @@ -1,63 +1,88 @@ Tests: - # Test 1: - - name : "Test 1: Simple mesh with flat topography" + + - name : "Test 1: Simple mesh with flat topography (P_SV wave)" description: > Testing mesh reader on a simple mesh with flat topography. config: nproc : 1 + elastic_wave: "P_SV" databases: mesh : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/database.bin" sources : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/sources.yaml" stations : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/STATIONS" - # Test 2: - - name : "Test 2: Simple mesh with curved topography" + + - name : "Test 2: Simple mesh with flat topography (SH wave)" + description: > + Testing mesh reader on a simple mesh with flat topography. + config: + nproc : 1 + elastic_wave: "SH" + databases: + mesh : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/database.bin" + sources : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/sources.yaml" + stations : "../../../tests/unit-tests/data/mesh/simple_mesh_flat_topography/STATIONS" + + - name : "Test 3: Simple mesh with curved topography" description: > Testing mesh reader on a simple mesh with curved topography. config: nproc : 1 + elastic_wave: "P_SV" databases: mesh : "../../../tests/unit-tests/data/mesh/simple_mesh_curved_topography/database.bin" sources : "../../../tests/unit-tests/data/mesh/simple_mesh_curved_topography/sources.yaml" stations : "../../../tests/unit-tests/data/mesh/simple_mesh_curved_topography/STATIONS" - # Test 3: - - name : "Test 3: Simple mesh with flat ocean bottom" + + - name : "Test 4: Simple mesh with flat ocean bottom" description: > Testing mesh reader on a simple mesh with flat ocean bottom. Fluid-solid interface is at the ocean bottom. config: nproc : 1 + elastic_wave: "P_SV" databases: mesh : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_flat_ocean_bottom/database.bin" sources : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_flat_ocean_bottom/sources.yaml" stations : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_flat_ocean_bottom/STATIONS" - # Test 4: - - name : "Test 4: Simple mesh with curved ocean bottom" + + - name : "Test 5: Simple mesh with curved ocean bottom" description: > Testing mesh reader on a simple mesh with curved ocean bottom. Fluid-solid interface is at the ocean bottom. config: nproc : 1 + elastic_wave: "P_SV" databases: mesh : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_curved_ocean_bottom/database.bin" sources : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_curved_ocean_bottom/sources.yaml" stations : "../../../tests/unit-tests/data/mesh/fluid_solid_mesh_curved_ocean_bottom/STATIONS" - # Test 5: - - name : "Test 5: Gmesh Example" + - name : "Test 6: Gmesh Example" description: > Testing mesh reader on a Gmesh example mesh. config: nproc : 1 + elastic_wave: "P_SV" databases: mesh : "../../../tests/unit-tests/data/mesh/Gmesh_Example_Stacey/database.bin" sources : "../../../tests/unit-tests/data/mesh/Gmesh_Example_Stacey/sources.yaml" stations : "../../../tests/unit-tests/data/mesh/Gmesh_Example_Stacey/STATIONS" + - name : "Test 7: Homogeneous Elastic Anisotropic Material (P_SV wave)" + description: > + Testing the reading of a homogeneous elastic anisotropic material. + config: + nproc : 1 + elastic_wave: "P_SV" + databases: + mesh : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/database.bin" + sources : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml" + stations : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/STATIONS" -# Test 6: - - name : "Test 6: Homogeneous Elastic Anisotropic Material" + - name : "Test 8: Homogeneous Elastic Anisotropic Material (SH wave)" description: > Testing the reading of a homogeneous elastic anisotropic material. config: nproc : 1 + elastic_wave: "SH" databases: mesh : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/database.bin" sources : "../../../tests/unit-tests/data/mesh/homogeneous_elastic_anisotropic/sources.yaml" diff --git a/tests/unit-tests/mesh/test_fixture/test_fixture.cpp b/tests/unit-tests/mesh/test_fixture/test_fixture.cpp index f5b9e8453..670299517 100644 --- a/tests/unit-tests/mesh/test_fixture/test_fixture.cpp +++ b/tests/unit-tests/mesh/test_fixture/test_fixture.cpp @@ -26,7 +26,9 @@ MESH::MESH() { for (auto &Test : Tests) { const auto [database_file, sources_file, stations_file] = Test.get_databases(); - specfem::mesh::mesh mesh = specfem::IO::read_mesh(database_file, mpi); + + const auto wave = Test.get_elastic_wave(); + specfem::mesh::mesh mesh = specfem::IO::read_mesh(database_file, wave, mpi); meshes.push_back(mesh); } diff --git a/tests/unit-tests/mesh/test_fixture/test_fixture.hpp b/tests/unit-tests/mesh/test_fixture/test_fixture.hpp index 9025ad713..f57de1b5a 100644 --- a/tests/unit-tests/mesh/test_fixture/test_fixture.hpp +++ b/tests/unit-tests/mesh/test_fixture/test_fixture.hpp @@ -32,18 +32,48 @@ struct database { std::string stations; }; +struct config { +public: + config() : nproc(1), elastic_wave("P_SV"){}; + config(const YAML::Node &Node) { + nproc = Node["nproc"].as(); + elastic_wave = Node["elastic_wave"].as(); + } + + int get_nproc() { return nproc; } + specfem::enums::elastic_wave get_elastic_wave() { + if (elastic_wave == "P_SV") + return specfem::enums::elastic_wave::p_sv; + else if (elastic_wave == "SH") + return specfem::enums::elastic_wave::sh; + else + throw std::runtime_error("Elastic wave type not supported"); + } + +private: + int nproc; + std::string elastic_wave; +}; + struct Test { public: Test(const YAML::Node &Node) { name = Node["name"].as(); description = Node["description"].as(); YAML::Node databases = Node["databases"]; + YAML::Node configuration = Node["config"]; try { database = test_configuration::database(databases); } catch (std::runtime_error &e) { throw std::runtime_error("Error in test configuration: " + name + "\n" + e.what()); } + try { + config = test_configuration::config(configuration); + } catch (std::runtime_error &e) { + throw std::runtime_error("Error in test configuration: " + name + "\n" + + e.what()); + } return; } @@ -51,9 +81,16 @@ struct Test { return database.get_databases(); } + int get_nproc() { return config.get_nproc(); } + + specfem::enums::elastic_wave get_elastic_wave() { + return config.get_elastic_wave(); + } + std::string name; std::string description; test_configuration::database database; + test_configuration::config config; }; } // namespace test_configuration