From ad4bcb6f1f0dc3c2f1340bb955f3bfee23b8795a Mon Sep 17 00:00:00 2001 From: Lucas Sawade Date: Wed, 26 Feb 2025 17:03:16 -0500 Subject: [PATCH] Updated fortran save_databases to write integer arrays and element_types to read integer arrays plus the first update for printing correcltuy --- .../generate_databases/save_arrays_solver.F90 | 58 +++++++++--- .../IO/mesh/impl/fortran/dim3/interface.hpp | 12 +++ .../mesh/dim3/element_types/element_types.hpp | 39 +++++--- src/IO/mesh/impl/fortran/dim3/mesh.cpp | 49 +++++----- src/mesh/dim3/element_types/element_types.cpp | 94 +++++++++++-------- 5 files changed, 163 insertions(+), 89 deletions(-) diff --git a/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 b/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 index 22bc4471..9a334ad8 100644 --- a/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 +++ b/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 @@ -110,30 +110,49 @@ end subroutine save_boundary_arrays ! subroutine that converts a logical array to an integer array and writes it ! a file IOUT (see save_global_arrays) - subroutine save_logical_nspec_array(array) + subroutine save_ispec_is_arrays(nspec, ispec_is_acoustic, ispec_is_elastic, ispec_is_poroelastic) use constants, only: IOUT - use generate_databases_par, only: nspec => NSPEC_AB - logical, dimension(nspec), intent(in) :: array + implicit none + + integer, intent(in) :: nspec + logical, dimension(nspec), intent(in) :: ispec_is_acoustic, ispec_is_elastic, ispec_is_poroelastic + ! To be written to file integer, dimension(nspec) :: array_int integer :: i + ! To be written to file + integer :: nspec_acoustic = 0, nspec_elastic = 0, nspec_poroelastic = 0 + + ! Loop over the number of element types in the array and convert the logical + ! values to integer values do i = 1, nspec - if (array(i)) then + array_int(i) = 0 + if (ispec_is_acoustic(i)) then + array_int(i) = 0 + nspec_acoustic = nspec_acoustic + 1 + elseif (ispec_is_elastic(i)) then array_int(i) = 1 + nspec_elastic = nspec_elastic + 1 + elseif (ispec_is_poroelastic(i)) then + array_int(i) = 2 + nspec_poroelastic = nspec_poroelastic + 1 else - array_int(i) = 0 - endif - enddo - + write(*,*) 'Error: ispec_is arrays are not correct' + stop + end if + end do + write(IOUT) nspec_acoustic + write(IOUT) nspec_elastic + write(IOUT) nspec_poroelastic write(IOUT) array_int - end subroutine save_logical_nspec_array + end subroutine save_ispec_is_arrays - end module save_arrays_module +end module save_arrays_module ! for external mesh subroutine save_arrays_solver_mesh() @@ -142,7 +161,7 @@ subroutine save_arrays_solver_mesh() use save_arrays_module, only: & save_global_arrays, & save_global_arrays_with_components, & - save_logical_nspec_array, & + save_ispec_is_arrays, & save_boundary_arrays use shared_parameters, only: ACOUSTIC_SIMULATION, ELASTIC_SIMULATION, POROELASTIC_SIMULATION, & @@ -274,10 +293,19 @@ subroutine save_arrays_solver_mesh() call save_global_arrays(nspec, kappastore) call save_global_arrays(nspec, mustore) - ! Converting the logical arrays to integer arrays and writing them as ints - write(IOUT) ispec_is_acoustic - write(IOUT) ispec_is_elastic - write(IOUT) ispec_is_poroelastic + ! save_ispec_is_arrays writes the following: + ! the number of each type of element + ! nspec_acoustic + ! nspec_elastic + ! nspec_poroelastic + ! + ! single integer array of nspec + ! represents the material type of each element + ! 0 = acoustic, 1 = elastic, 2 = poroelastic + + call save_ispec_is_arrays(nspec, ispec_is_acoustic, & + ispec_is_elastic, & + ispec_is_poroelastic) ! stamp for checking i/o itest = 9999 diff --git a/include/IO/mesh/impl/fortran/dim3/interface.hpp b/include/IO/mesh/impl/fortran/dim3/interface.hpp index 60e95ef0..f9e77ba1 100644 --- a/include/IO/mesh/impl/fortran/dim3/interface.hpp +++ b/include/IO/mesh/impl/fortran/dim3/interface.hpp @@ -33,6 +33,18 @@ void read_ibool(std::ifstream &stream, specfem::mesh::mapping &mapping, const specfem::MPI::MPI *mpi); +/** + * @brief Read element types from 3D mesh database + * + * @param stream Input stream + * @param element_types Element types object + * @param mpi MPI object + */ +void read_element_types( + std::ifstream &stream, + specfem::mesh::element_types &element_types, + const specfem::MPI::MPI *mpi); + /** * @brief Read coordinates from 3D mesh database * diff --git a/include/mesh/dim3/element_types/element_types.hpp b/include/mesh/dim3/element_types/element_types.hpp index 3c5cf71f..3594eddd 100644 --- a/include/mesh/dim3/element_types/element_types.hpp +++ b/include/mesh/dim3/element_types/element_types.hpp @@ -19,8 +19,8 @@ template <> struct element_types { int nspec; ///< Number of spectral elements - int nelastic; ///< Number of elastic spectral elements int nacoustic; ///< Number of acoustic spectral elements + int nelastic; ///< Number of elastic spectral elements int nporoelastic; ///< Number of poroelastic spectral elements template using View1D = Kokkos::View; @@ -40,6 +40,9 @@ template <> struct element_types { * @brief Construct a new element types object * * @param nspec Number of spectral elements + * @param nacoustic Number of acoustic spectral elements + * @param nelastic Number of elastic spectral elements + * @param nporoelastic Number of poroelastic spectral elements * * @code{.cpp} * // Example of how to use this constructor @@ -48,10 +51,14 @@ template <> struct element_types { * element_types(nspec); * @endcode */ - element_types(const int nspec) - : nspec(nspec), ispec_is_elastic("ispec_is_elastic", nspec), - ispec_is_acoustic("ispec_is_acoustic", nspec), - ispec_is_poroelastic("ispec_is_poroelastic", nspec){}; + element_types(const int nspec, const int nacoustic, const int nelastic, + const int nporoelastic) + : nspec(nspec), nacoustic(nacoustic), nelastic(nelastic), + nporoelastic(nporoelastic), ispec_type("ispec_type", nspec), + ispec_acoustic("mesh.element_types.ispec_acoustic", nacoustic), + ispec_elastic("mesh.element_types.ispec_elastic", nelastic), + ispec_poroelastic("mesh.element_types.ispec_poroelastic", + nporoelastic){}; ///@} @@ -65,7 +72,7 @@ template <> struct element_types { * * @see get_elements */ - void set_elements(); + void set_elements(View1D &ispec_type_in); /** * @brief Get the elements for the given medium type @@ -79,14 +86,14 @@ template <> struct element_types { * @brief Print basic information about the element types * */ - void print() const; + std::string print() const; /** * @brief Print the element type of a specific spectral element * * @param ispec Index of the spectral element */ - void print(const int ispec) const; + std::string print(const int ispec) const; /** * @brief Print the element type of element i/n @@ -94,16 +101,18 @@ template <> struct element_types { * @tparam MediumTag Medium type */ template - void print(const int i) const; + std::string print(const int i) const; - View1D ispec_is_elastic; ///< Elastic spectral elements - View1D ispec_is_acoustic; ///< Acoustic spectral elements - View1D ispec_is_poroelastic; ///< Poroelastic spectral elements + View1D + ispec_type; ///< Elastic spectral elements with: + ///< ispec_type[ispec] = 0 (acoustic) + ///< ispec_type[ispec] = 1 (elastic) + ///< ispec_type[ispec] = 2 (poroelastic) private: - View1D ispec_elastic; ///< Elastic spectral elements - View1D ispec_acoustic; ///< Acoustic spectral elements - View1D ispec_poroelastic; ///< Poroelastic spectral elements + View1D ispec_elastic; ///< Elastic to global mapping + View1D ispec_acoustic; ///< Acoustic to global mapping + View1D ispec_poroelastic; ///< Poroelastic to global mapping }; } // namespace mesh } // namespace specfem diff --git a/src/IO/mesh/impl/fortran/dim3/mesh.cpp b/src/IO/mesh/impl/fortran/dim3/mesh.cpp index 8ba53bda..422357dd 100644 --- a/src/IO/mesh/impl/fortran/dim3/mesh.cpp +++ b/src/IO/mesh/impl/fortran/dim3/mesh.cpp @@ -22,18 +22,21 @@ #include template -void try_print_medium_element( +std::string try_print_medium_element( const specfem::mesh::element_types &elements_types, int index) { + + std::ostringstream message; try { - elements_types.print(index); + message << elements_types.print(index); } catch (std::runtime_error &e) { - std::cout << e.what(); + message << e.what(); } catch (...) { - std::cout << "Unknown exception caught in try_print_medium_element" - << std::endl; + message << "Unknown exception caught in try_print_medium_element" + << ".\n"; } + return message.str(); } specfem::mesh::mesh @@ -197,34 +200,36 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, mesh.materials.print(); #endif + int nacoustic, nelastic, nporoelastic; + + try_read_line("read_nacoustic", stream, &nacoustic); + try_read_line("read_nelastic", stream, &nelastic); + try_read_line("read_nporoelastic", stream, &nporoelastic); + // Initialize element types object mesh.elements_types = specfem::mesh::element_types( - mesh.parameters.nspec); + mesh.parameters.nspec, nacoustic, nelastic, nporoelastic); - // Read boolean is x array! - try_read_array("read_ispec_is_acoustic", stream, - mesh.elements_types.ispec_is_acoustic); - try_read_array("read_ispec_is_elastic", stream, - mesh.elements_types.ispec_is_elastic); - try_read_array("read_ispec_is_poroelastic", stream, - mesh.elements_types.ispec_is_poroelastic); + // Reading the ispec type + auto ispec_type = Kokkos::View("ispec_type", nspec); + try_read_array("read_ispec_type", stream, ispec_type); // Compute the element arrays - mesh.elements_types.set_elements(); + mesh.elements_types.set_elements(ispec_type); #ifndef NDEBUG // Print the element types - mesh.elements_types.print(); - mesh.elements_types.print(0); + mpi->cout(mesh.elements_types.print()); + mpi->cout(mesh.elements_types.print(0)); // Print elements first element of each category - try_print_medium_element( - mesh.elements_types, 0); - try_print_medium_element( - mesh.elements_types, 0); - try_print_medium_element( - mesh.elements_types, 0); + mpi->cout(try_print_medium_element( + mesh.elements_types, 0)); + mpi->cout(try_print_medium_element( + mesh.elements_types, 0)); + mpi->cout(try_print_medium_element( + mesh.elements_types, 0)); #endif // Read test value 9999 diff --git a/src/mesh/dim3/element_types/element_types.cpp b/src/mesh/dim3/element_types/element_types.cpp index d996e732..f75c70cb 100644 --- a/src/mesh/dim3/element_types/element_types.cpp +++ b/src/mesh/dim3/element_types/element_types.cpp @@ -1,14 +1,15 @@ #include "mesh/dim3/element_types/element_types.hpp" #include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" #include -void specfem::mesh::element_types< - specfem::dimension::type::dim3>::set_elements() { +void specfem::mesh::element_types::set_elements( + Kokkos::View &ispec_type_in) { // Create the vectors - std::vector ispec_elastic; - std::vector ispec_acoustic; - std::vector ispec_poroelastic; + // std::vector ispec_elastic; + // std::vector ispec_acoustic; + // std::vector ispec_poroelastic; // Initialize the number of elements this->nelastic = 0; @@ -17,14 +18,17 @@ void specfem::mesh::element_types< // Initialize the vectors for (int ispec = 0; ispec < nspec; ispec++) { - if (ispec_is_elastic(ispec) == true) { - ispec_elastic.push_back(ispec); - this->nelastic++; - } else if (ispec_is_acoustic(ispec) == true) { - ispec_acoustic.push_back(ispec); + if (ispec_type_in[ispec] == 0) { + this->ispec_type(ispec) = specfem::element::medium_tag::acoustic; + this->ispec_acoustic(nacoustic) = ispec; this->nacoustic++; - } else if (ispec_is_poroelastic(ispec) == true) { - ispec_poroelastic.push_back(ispec); + } else if (ispec_type_in[ispec] == 1) { + this->ispec_type(ispec) = specfem::element::medium_tag::elastic; + this->ispec_elastic(nelastic) = ispec; + this->nelastic++; + } else if (ispec_type_in[ispec] == 2) { + this->ispec_type(ispec) = specfem::element::medium_tag::poroelastic; + this->ispec_poroelastic(nporoelastic) = ispec; this->nporoelastic++; } else { std::ostringstream message; @@ -32,35 +36,46 @@ void specfem::mesh::element_types< << " is not assigned to any material"; throw std::runtime_error(message.str()); } - - // Assign the vectors to the view - this->ispec_elastic = View1D("ispec_elastic", ispec_elastic.size()); - this->ispec_acoustic = View1D("ispec_acoustic", ispec_acoustic.size()); - this->ispec_poroelastic = - View1D("ispec_poroelastic", ispec_poroelastic.size()); }; } -void specfem::mesh::element_types::print() - const { - std::cout << "Number of elastic elements: " << nelastic << std::endl; - std::cout << "Number of acoustic elements: " << nacoustic << std::endl; - std::cout << "Number of poroelastic elements: " << nporoelastic << std::endl; +std::string +specfem::mesh::element_types::print() const { + + std::ostringstream message; + message << "Number of acoustic elements: " << nacoustic << ".\n"; + message << "Number of elastic elements: " << nelastic << ".\n"; + message << "Number of poroelastic elements: " << nporoelastic << ".\n"; + + return message.str(); } -void specfem::mesh::element_types::print( +std::string specfem::mesh::element_types::print( const int ispec) const { - std::cout << "Element " << ispec << " is elastic: " << ispec_is_elastic(ispec) - << std::endl; - std::cout << "Element " << ispec - << " is acoustic: " << ispec_is_acoustic(ispec) << std::endl; - std::cout << "Element " << ispec - << " is poroelastic: " << ispec_is_poroelastic(ispec) << std::endl; + + std::ostringstream message; + + if (ispec >= nspec) { + std::ostringstream message; + message << "Error: Spectral element " << ispec + << " does not exist.\nNumber of spectral elements is :" << nspec + << ". " + << "(" << __FILE__ << ":" << __LINE__ << ")\n"; + throw std::runtime_error(message.str()); + } + + message << "Element " << ispec << " is "; + message << specfem::element::to_string(ispec_type(ispec)) << ".\n"; + + return message.str(); } template -void specfem::mesh::element_types::print( +std::string specfem::mesh::element_types::print( const int i) const { + + std::ostringstream message; + if (MediumTag == specfem::element::medium_tag::elastic) { if (i >= nelastic) { std::ostringstream message; @@ -71,7 +86,8 @@ void specfem::mesh::element_types::print( << "(" << __FILE__ << ":" << __LINE__ << ")\n"; throw std::runtime_error(message.str()); } - std::cout << "Elastic element " << ispec_elastic(i) << std::endl; + message << "Elastic element " << i << " is global element " + << ispec_elastic(i) << ".\n"; } else if (MediumTag == specfem::element::medium_tag::acoustic) { if (i >= nacoustic) { @@ -84,7 +100,8 @@ void specfem::mesh::element_types::print( throw std::runtime_error(message.str()); } - std::cout << "Acoustic element " << ispec_acoustic(i) << std::endl; + message << "Acoustic element " << i << " is global element " + << ispec_acoustic(i) << ".\n"; } else if (MediumTag == specfem::element::medium_tag::poroelastic) { @@ -98,17 +115,20 @@ void specfem::mesh::element_types::print( throw std::runtime_error(message.str()); } - std::cout << "Poroelastic element " << ispec_poroelastic(i) << std::endl; + message << "Poroelastic element " << i << " is global element " + << ispec_poroelastic(i) << ".\n"; } + + return message.str(); } // Explicit instantiations -template void +template std::string specfem::mesh::element_types::print< specfem::element::medium_tag::elastic>(const int) const; -template void +template std::string specfem::mesh::element_types::print< specfem::element::medium_tag::acoustic>(const int) const; -template void +template std::string specfem::mesh::element_types::print< specfem::element::medium_tag::poroelastic>(const int) const;