Skip to content

Commit

Permalink
Updated fortran save_databases to write integer arrays and element_ty…
Browse files Browse the repository at this point in the history 
…pes to read integer arrays plus the first update for printing correcltuy
  • Loading branch information
@lsawade
lsawade committed Feb 26, 2025
1 parent 3af4632 commit ad4bcb6
Show file tree
Hide file tree
Showing 5 changed files with 163 additions and 89 deletions.
58 changes: 43 additions & 15 deletions fortran/meshfem3d/generate_databases/save_arrays_solver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand All @@ -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, &
Expand Down Expand Up @@ -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
Expand Down
12 changes: 12 additions & 0 deletions include/IO/mesh/impl/fortran/dim3/interface.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,6 +33,18 @@ void read_ibool(std::ifstream &stream,
specfem::mesh::mapping<specfem::dimension::type::dim3> &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<specfem::dimension::type::dim3> &element_types,
const specfem::MPI::MPI *mpi);

/**
* @brief Read coordinates from 3D mesh database
*
Expand Down
39 changes: 24 additions & 15 deletions include/mesh/dim3/element_types/element_types.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ template <> struct element_types<specfem::dimension::type::dim3> {

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 <typename T> using View1D = Kokkos::View<T *, Kokkos::HostSpace>;
Expand All @@ -40,6 +40,9 @@ template <> struct element_types<specfem::dimension::type::dim3> {
* @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
Expand All @@ -48,10 +51,14 @@ template <> struct element_types<specfem::dimension::type::dim3> {
* 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){};

///@}

Expand All @@ -65,7 +72,7 @@ template <> struct element_types<specfem::dimension::type::dim3> {
*
* @see get_elements
*/
void set_elements();
void set_elements(View1D<int> &ispec_type_in);

/**
* @brief Get the elements for the given medium type
Expand All @@ -79,31 +86,33 @@ template <> struct element_types<specfem::dimension::type::dim3> {
* @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<medium>
*
* @tparam MediumTag Medium type
*/
template <specfem::element::medium_tag MediumTag>
void print(const int i) const;
std::string print(const int i) const;

View1D<bool> ispec_is_elastic; ///< Elastic spectral elements
View1D<bool> ispec_is_acoustic; ///< Acoustic spectral elements
View1D<bool> ispec_is_poroelastic; ///< Poroelastic spectral elements
View1D<specfem::element::medium_tag>
ispec_type; ///< Elastic spectral elements with:
///< ispec_type[ispec] = 0 (acoustic)
///< ispec_type[ispec] = 1 (elastic)
///< ispec_type[ispec] = 2 (poroelastic)

private:
View1D<int> ispec_elastic; ///< Elastic spectral elements
View1D<int> ispec_acoustic; ///< Acoustic spectral elements
View1D<int> ispec_poroelastic; ///< Poroelastic spectral elements
View1D<int> ispec_elastic; ///< Elastic to global mapping
View1D<int> ispec_acoustic; ///< Acoustic to global mapping
View1D<int> ispec_poroelastic; ///< Poroelastic to global mapping
};
} // namespace mesh
} // namespace specfem
49 changes: 27 additions & 22 deletions src/IO/mesh/impl/fortran/dim3/mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,18 +22,21 @@
#include <vector>

template <specfem::element::medium_tag medium>
void try_print_medium_element(
std::string try_print_medium_element(
const specfem::mesh::element_types<specfem::dimension::type::dim3>
&elements_types,
int index) {

std::ostringstream message;
try {
elements_types.print<medium>(index);
message << elements_types.print<medium>(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<specfem::dimension::type::dim3>
Expand Down Expand Up @@ -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<specfem::dimension::type::dim3>(
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<int *, Kokkos::HostSpace>("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<specfem::element::medium_tag::acoustic>(
mesh.elements_types, 0);
try_print_medium_element<specfem::element::medium_tag::elastic>(
mesh.elements_types, 0);
try_print_medium_element<specfem::element::medium_tag::poroelastic>(
mesh.elements_types, 0);
mpi->cout(try_print_medium_element<specfem::element::medium_tag::acoustic>(
mesh.elements_types, 0));
mpi->cout(try_print_medium_element<specfem::element::medium_tag::elastic>(
mesh.elements_types, 0));
mpi->cout(try_print_medium_element<specfem::element::medium_tag::poroelastic>(
mesh.elements_types, 0));
#endif

// Read test value 9999
Expand Down
Loading

0 comments on commit ad4bcb6

Please sign in to comment.