diff --git a/CMakeLists.txt b/CMakeLists.txt index 5c7fa6ee4..3510603ca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -245,6 +245,7 @@ add_library( src/IO/mesh/impl/fortran/dim3/read_parameters.cpp src/IO/mesh/impl/fortran/dim3/read_coordinates.cpp src/IO/mesh/impl/fortran/dim3/read_partial_derivatives.cpp + src/IO/mesh/impl/fortran/dim3/utilities.cpp ) if (NOT HDF5_CXX_BUILD) @@ -350,15 +351,22 @@ add_library( src/mesh/dim2/tags/tags.cpp # 3-D src/mesh/dim3/mesh.cpp + src/mesh/dim3/boundaries/absorbing_boundary.cpp + src/mesh/dim3/boundaries/free_surface.cpp + src/mesh/dim3/coupled_interfaces/coupled_interfaces.cpp + src/mesh/dim3/element_types/element_types.cpp src/mesh/dim3/parameters/parameters.cpp src/mesh/dim3/parameters/parameters.cpp src/mesh/dim3/mapping/mapping.cpp + src/mesh/dim3/materials/materials.cpp src/mesh/dim3/coordinates/coordinates.cpp src/mesh/dim3/partial_derivatives/partial_derivatives.cpp + ) target_link_libraries( mesh + enumerations Kokkos::kokkos specfem_mpi # material_class diff --git a/docs/Doxyfile.in b/docs/Doxyfile.in index 08946ddfc..975b10659 100644 --- a/docs/Doxyfile.in +++ b/docs/Doxyfile.in @@ -502,7 +502,7 @@ NUM_PROC_THREADS = 1 # normally produced when WARNINGS is set to YES. # The default value is: NO. -EXTRACT_ALL = NO +EXTRACT_ALL = YES # If the EXTRACT_PRIVATE tag is set to YES, all private members of a class will # be included in the documentation. @@ -741,7 +741,7 @@ GENERATE_DEPRECATEDLIST= YES # sections, marked by \if ... \endif and \cond # ... \endcond blocks. -ENABLED_SECTIONS = +ENABLED_SECTIONS = NO # The MAX_INITIALIZER_LINES tag determines the maximum number of lines that the # initial value of a variable or macro / define can have for it to appear in the @@ -1222,7 +1222,7 @@ IGNORE_PREFIX = # If the GENERATE_HTML tag is set to YES, doxygen will generate HTML output # The default value is: YES. -GENERATE_HTML = YES +GENERATE_HTML = NO # The HTML_OUTPUT tag is used to specify where the HTML docs will be put. If a # relative path is entered the value of OUTPUT_DIRECTORY will be put in front of @@ -1230,7 +1230,7 @@ GENERATE_HTML = YES # The default directory is: html. # This tag requires that the tag GENERATE_HTML is set to YES. -HTML_OUTPUT = html +HTML_OUTPUT = "doxygen/html_doxygen" # The HTML_FILE_EXTENSION tag can be used to specify the file extension for each # generated HTML page (for example: .htm, .php, .asp). @@ -2146,7 +2146,7 @@ GENERATE_XML = YES # The default directory is: xml. # This tag requires that the tag GENERATE_XML is set to YES. -XML_OUTPUT = xml +XML_OUTPUT = doxygen/xml # If the XML_PROGRAMLISTING tag is set to YES, doxygen will dump the program # listings (including syntax highlighting and cross-referencing information) to diff --git a/docs/api/IO/fortran_io.rst b/docs/api/IO/fortran_io.rst index 86982c3b8..47e590ea8 100644 --- a/docs/api/IO/fortran_io.rst +++ b/docs/api/IO/fortran_io.rst @@ -5,4 +5,4 @@ Fortran IO Fortran IO namespace provides a functions to read data stored in unformatted fortran binary files. -.. doxygenfunction:: specfem::fortran_IO::fortran_read_line +.. doxygenfunction:: specfem::IO::fortran_read_line diff --git a/docs/api/IO/index.rst b/docs/api/IO/index.rst index 766c15f66..407a17d6d 100644 --- a/docs/api/IO/index.rst +++ b/docs/api/IO/index.rst @@ -29,8 +29,11 @@ are exposed through the following modules: reader/index -Read Mesh ---------- +Input for the 2D Solver ++++++++++++++++++++++++ + +Read the 2D Mesh +---------------- .. doxygenfunction:: specfem::IO::read_2d_mesh @@ -45,3 +48,35 @@ Read Receivers -------------- .. doxygenfunction:: specfem::IO::read_receivers(const std::string stations_file, const type_real angle) + + + + +Input for the 3D Solver ++++++++++++++++++++++++ + + +Read the 3D Mesh +---------------- + +.. doxygenfunction:: specfem::IO::read_3d_mesh + + +Helper functions +'''''''''''''''' + +**Reading any values** + +.. doxygenfunction:: specfem::IO::mesh::impl::fortran::dim3::try_read_line + +**Reading any array** + +.. doxygenfunction:: specfem::IO::mesh::impl::fortran::dim3::try_read_array + +.. doxygenfunction:: specfem::IO::mesh::impl::fortran::dim3::read_array + +**Reading index arrays** + +.. doxygenfunction:: specfem::IO::mesh::impl::fortran::dim3::try_read_index_array + +.. doxygenfunction:: specfem::IO::mesh::impl::fortran::dim3::read_index_array diff --git a/docs/api/enumerations/index.rst b/docs/api/enumerations/index.rst index fb0d34186..d7b7184eb 100644 --- a/docs/api/enumerations/index.rst +++ b/docs/api/enumerations/index.rst @@ -1,4 +1,4 @@ -.. _enumerations:: +.. _enumerations: Enumerations ============ diff --git a/docs/api/medium/index.rst b/docs/api/medium/index.rst index c90872a8c..e1d3c0797 100644 --- a/docs/api/medium/index.rst +++ b/docs/api/medium/index.rst @@ -1,4 +1,4 @@ -.. _medium:: +.. _medium: Medium modules ============== diff --git a/docs/api/mesh/boundaries/absorbing_boundaries.rst b/docs/api/mesh/dim2/boundaries/absorbing_boundaries.rst similarity index 100% rename from docs/api/mesh/boundaries/absorbing_boundaries.rst rename to docs/api/mesh/dim2/boundaries/absorbing_boundaries.rst diff --git a/docs/api/mesh/boundaries/acoustic_free_surface.rst b/docs/api/mesh/dim2/boundaries/acoustic_free_surface.rst similarity index 100% rename from docs/api/mesh/boundaries/acoustic_free_surface.rst rename to docs/api/mesh/dim2/boundaries/acoustic_free_surface.rst diff --git a/docs/api/mesh/boundaries/index.rst b/docs/api/mesh/dim2/boundaries/index.rst similarity index 100% rename from docs/api/mesh/boundaries/index.rst rename to docs/api/mesh/dim2/boundaries/index.rst diff --git a/docs/api/mesh/control_nodes/index.rst b/docs/api/mesh/dim2/control_nodes/index.rst similarity index 100% rename from docs/api/mesh/control_nodes/index.rst rename to docs/api/mesh/dim2/control_nodes/index.rst diff --git a/docs/api/mesh/coupled_interfaces/index.rst b/docs/api/mesh/dim2/coupled_interfaces/index.rst similarity index 69% rename from docs/api/mesh/coupled_interfaces/index.rst rename to docs/api/mesh/dim2/coupled_interfaces/index.rst index 92194b71d..15413cfe1 100644 --- a/docs/api/mesh/coupled_interfaces/index.rst +++ b/docs/api/mesh/dim2/coupled_interfaces/index.rst @@ -4,7 +4,7 @@ Coupled Interfaces ================== -.. doxygenstruct:: specfem::mesh::coupled_interfaces +.. doxygenstruct:: specfem::mesh::coupled_interfaces< specfem::dimension::type::dim2 > :members: diff --git a/docs/api/mesh/dim2/index.rst b/docs/api/mesh/dim2/index.rst new file mode 100644 index 000000000..2f4b8efd4 --- /dev/null +++ b/docs/api/mesh/dim2/index.rst @@ -0,0 +1,18 @@ + +.. _mesh_2D: + +SPECFEM++ 2D Mesh +================= + +.. doxygenstruct:: specfem::mesh::mesh< specfem::dimension::type::dim2 > + :members: + + +.. toctree:: + :maxdepth: 1 + + boundaries/index + control_nodes/index + coupled_interfaces/index + materials/index + tags/index diff --git a/docs/api/mesh/materials/index.rst b/docs/api/mesh/dim2/materials/index.rst similarity index 100% rename from docs/api/mesh/materials/index.rst rename to docs/api/mesh/dim2/materials/index.rst diff --git a/docs/api/mesh/tags/index.rst b/docs/api/mesh/dim2/tags/index.rst similarity index 100% rename from docs/api/mesh/tags/index.rst rename to docs/api/mesh/dim2/tags/index.rst diff --git a/docs/api/mesh/dim3/boundaries/index.rst b/docs/api/mesh/dim3/boundaries/index.rst new file mode 100644 index 000000000..ae1efa801 --- /dev/null +++ b/docs/api/mesh/dim3/boundaries/index.rst @@ -0,0 +1,23 @@ + +.. _mesh_boundaries_3D: + +Boundaries +++++++++++ + + +.. doxygenstruct:: specfem::mesh::boundaries< specfem::dimension::type::dim3 > + :members: + +Absorbing Boundaries +==================== + +.. doxygenstruct:: specfem::mesh::absorbing_boundary< specfem::dimension::type::dim3 > + :members: + + + +Free surface boundary +===================== + +.. doxygenstruct:: specfem::mesh::free_surface< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/coordinates/index.rst b/docs/api/mesh/dim3/coordinates/index.rst new file mode 100644 index 000000000..691b1a4f5 --- /dev/null +++ b/docs/api/mesh/dim3/coordinates/index.rst @@ -0,0 +1,7 @@ +.. _mesh_coordinates_3D: + +Coordinates ++++++++++++ + +.. doxygenstruct:: specfem::mesh::coordinates< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/coupled_interfaces/index.rst b/docs/api/mesh/dim3/coupled_interfaces/index.rst new file mode 100644 index 000000000..472388c1b --- /dev/null +++ b/docs/api/mesh/dim3/coupled_interfaces/index.rst @@ -0,0 +1,8 @@ +.. mesh_3D_coupled_interfaces: + +Coupled Interfaces +================== + + +.. doxygenstruct:: specfem::mesh::coupled_interfaces< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/element_types/index.rst b/docs/api/mesh/dim3/element_types/index.rst new file mode 100644 index 000000000..ae87c7550 --- /dev/null +++ b/docs/api/mesh/dim3/element_types/index.rst @@ -0,0 +1,7 @@ +.. mesh_3D_element_types: + +Element Types +============= + +.. doxygenstruct:: specfem::mesh::element_types< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/index.rst b/docs/api/mesh/dim3/index.rst new file mode 100644 index 000000000..96cfc2154 --- /dev/null +++ b/docs/api/mesh/dim3/index.rst @@ -0,0 +1,25 @@ + +.. _mesh_3D: + +SPECFEM++ 3D Mesh +================= + +.. doxygenstruct:: specfem::mesh::mesh< specfem::dimension::type::dim3 > + :members: + + +What the 3D mesh struct holds +----------------------------- + +.. toctree:: + :maxdepth: 1 + + parameters/index + coordinates/index + mapping/index + partial_derivatives/index + element_types/index + materials/index + mass_matrix/index + boundaries/index + coupled_interfaces/index diff --git a/docs/api/mesh/dim3/mapping/index.rst b/docs/api/mesh/dim3/mapping/index.rst new file mode 100644 index 000000000..2311fe412 --- /dev/null +++ b/docs/api/mesh/dim3/mapping/index.rst @@ -0,0 +1,7 @@ +.. mesh_3D_mapping: + +Mapping +======= + +.. doxygenstruct:: specfem::mesh::mapping< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/mass_matrix/index.rst b/docs/api/mesh/dim3/mass_matrix/index.rst new file mode 100644 index 000000000..8bfe72024 --- /dev/null +++ b/docs/api/mesh/dim3/mass_matrix/index.rst @@ -0,0 +1,7 @@ +.. mesh_3D_mass_matrix: + +Matrix +====== + +.. doxygenstruct:: specfem::mesh::mass_matrix< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/materials/index.rst b/docs/api/mesh/dim3/materials/index.rst new file mode 100644 index 000000000..69d32d5d7 --- /dev/null +++ b/docs/api/mesh/dim3/materials/index.rst @@ -0,0 +1,7 @@ +.. mesh_3D_materials: + +Material +======== + +.. doxygenstruct:: specfem::mesh::materials< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/parameters/index.rst b/docs/api/mesh/dim3/parameters/index.rst new file mode 100644 index 000000000..18c727fd6 --- /dev/null +++ b/docs/api/mesh/dim3/parameters/index.rst @@ -0,0 +1,7 @@ +.. mesh_3D_parameters: + +Parameters +========== + +.. doxygenstruct:: specfem::mesh::parameters< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/dim3/partial_derivatives/index.rst b/docs/api/mesh/dim3/partial_derivatives/index.rst new file mode 100644 index 000000000..46adbc49d --- /dev/null +++ b/docs/api/mesh/dim3/partial_derivatives/index.rst @@ -0,0 +1,7 @@ +.. mesh_3D_partial_derivatives: + +Partial Derivatives +=================== + +.. doxygenstruct:: specfem::mesh::partial_derivatives< specfem::dimension::type::dim3 > + :members: diff --git a/docs/api/mesh/index.rst b/docs/api/mesh/index.rst index 8d864e45c..f67b2ecc8 100644 --- a/docs/api/mesh/index.rst +++ b/docs/api/mesh/index.rst @@ -11,8 +11,5 @@ SPECFEM++ Mesh .. toctree:: :maxdepth: 1 - boundaries/index - control_nodes/index - coupled_interfaces/index - materials/index - tags/index + dim2/index + dim3/index diff --git a/docs/conf.py b/docs/conf.py index 38557d484..ed14d1f71 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -14,9 +14,31 @@ # import sys # sys.path.insert(0, os.path.abspath('.')) import subprocess +import sys +import os # Doxygen -subprocess.call("doxygen Doxyfile.in", shell=True) +doxygen_cmd = "doxygen Doxyfile.in" # or Doxyfile.in if that's the correct filename + +# Create the doxygen subdirectory if it doesn't exist +if os.path.exists("_build/doxygen") is False: + os.makedirs("_build/doxygen") + +if os.environ.get("NODOXYGEN") is not None: + print("Skipping Doxygen build as NODOXYGEN environment variable is set") + result = 0 +else: + result = subprocess.call(doxygen_cmd, shell=True) + +if result != 0: + print(f"Error: Doxygen command '{doxygen_cmd}' failed with code {result}") + sys.exit(1) + +# Check if output directory exists +if not os.path.exists("_build/doxygen/xml"): + print("Error: Doxygen XML output directory '_build/doxygen/xml' not found") + print("Please check Doxygen configuration and make sure it runs correctly") + sys.exit(1) # -- Project information ----------------------------------------------------- @@ -42,10 +64,14 @@ "sphinx.ext.viewcode", "sphinx_sitemap", "sphinx.ext.inheritance_diagram", + "sphinx_design", "breathe", "sphinx_copybutton", ] +# Adding this to avoid the WARNING: duplicate label warning +autosectionlabel_prefix_document = True + # Add any paths that contain templates here, relative to this directory. templates_path = ["_templates"] @@ -78,7 +104,7 @@ # } # html_logo = '' github_url = "https://github.com/PrincetonUniversity/SPECFEMPP" -# html_baseurl = '' +html_baseurl = "https://specfem2d-kokkos.readthedocs.io/" # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, @@ -96,6 +122,6 @@ # -- Breathe configuration ------------------------------------------------- -breathe_projects = {"SPECFEM KOKKOS IMPLEMENTATION": "_build/xml"} +breathe_projects = {"SPECFEM KOKKOS IMPLEMENTATION": "_build/doxygen/xml"} breathe_default_project = "SPECFEM KOKKOS IMPLEMENTATION" breathe_default_members = () diff --git a/docs/developer_documentation/tutorials/tutorial1/Chapter1/index.rst b/docs/developer_documentation/tutorials/tutorial1/Chapter1/index.rst index e3236b6ee..88dcd5b22 100644 --- a/docs/developer_documentation/tutorials/tutorial1/Chapter1/index.rst +++ b/docs/developer_documentation/tutorials/tutorial1/Chapter1/index.rst @@ -172,7 +172,7 @@ Defining sources Now that we have the mesh and receiver locations, we can define the sources. We will use the following source file to define the sources. -.. code:: yaml +.. code-block:: yaml :caption: sources.yaml number-of-sources: 1 @@ -189,6 +189,8 @@ Now that we have the mesh and receiver locations, we can define the sources. We tshift: 0.0 f0: 10.0 + + Final mesh ---------- diff --git a/docs/meshfem2d/index.rst b/docs/meshfem2d/index.rst index 8b181cd7f..5212eea37 100644 --- a/docs/meshfem2d/index.rst +++ b/docs/meshfem2d/index.rst @@ -3,9 +3,18 @@ MESHFEM Parameter Documentation =============================== -SPECFEM++ uses a modified version of the mesher provided by the original `SPECFEM2D `_ code. The modifications isolate only the nacessary parameters for the meshing process and remove those needed by the solver from the original ``Par_File``. Here I will document the parameters that are used by the meshing process. However you should also refer to the `SPECFEM2D documentation `_ for a more in-depth description of the mesher. +SPECFEM++ uses a modified version of the mesher provided by the original +`SPECFEM2D `_ +code. The modifications isolate only the nacessary parameters for the meshing +process and remove those needed by the solver from the original ``Par_File``. +Here I will document the parameters that are used by the meshing process. +However you should also refer to the `SPECFEM2D documentation +`_ for a more +in-depth description of the mesher. -To define the meshing parameters, you will need to create a ``Par_File``. This is a simple text file that contains the parameters you wish to use. The parameters are defined in the following format: +To define the meshing parameters, you will need to create a ``Par_File``. This +is a simple text file that contains the parameters you wish to use. The +parameters are defined in the following format: .. code-block:: bash @@ -31,7 +40,11 @@ Parameter Description Topography File --------------- -Topography files are used to define the surface topography of the mesh. The topography file is a simple text file that describes the topography of every interface in the simulation domain. For example the following topography file describes a simple 2 layer model with a flat surface and a flat interface between the two layers: +Topography files are used to define the surface topography of the mesh. The +topography file is a simple text file that describes the topography of every +interface in the simulation domain. For example the following topography file +describes a simple 2 layer model with a flat surface and a flat interface +between the two layers: .. code-block:: bash diff --git a/docs/meshfem2d/meshing_parameters.rst b/docs/meshfem2d/meshing_parameters.rst index 8bf46669d..291792a14 100644 --- a/docs/meshfem2d/meshing_parameters.rst +++ b/docs/meshfem2d/meshing_parameters.rst @@ -3,7 +3,7 @@ Meshing Parameters ================== **Parameter Name**: ``PARTITIONING_TYPE`` -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. note:: @@ -17,7 +17,7 @@ Meshing Parameters **Type**: ``int`` **Parameter Name**: ``NGNOD`` -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **Description**: Number of control nodes per element. @@ -25,7 +25,7 @@ Meshing Parameters - Supported values: ``9`` **Parameter Name**: ``database_filename`` -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **Description**: Name of the database file to store the generated mesh. diff --git a/docs/meshfem2d/velocity_models.rst b/docs/meshfem2d/velocity_models.rst index 534e99214..8d1cb7ac0 100644 --- a/docs/meshfem2d/velocity_models.rst +++ b/docs/meshfem2d/velocity_models.rst @@ -3,7 +3,7 @@ Velocity Models ================ **Parameter Name**: ``nbmodels`` -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ **Description**: Number of material systems in the model. diff --git a/docs/parameter_documentation/databases.rst b/docs/parameter_documentation/databases.rst deleted file mode 100644 index 35ee27310..000000000 --- a/docs/parameter_documentation/databases.rst +++ /dev/null @@ -1,35 +0,0 @@ -Databases -######### - -Databases section defines location of database files. - -Parameter definitions -===================== - -**Parameter name** : ``databases`` -------------------------------- - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Define databases section - -.. _database-file-parameter: - -**Parameter name** : ``databases.mesh-database`` -****************************************************** - -**default value**: None - -**possible values**: [string] - -**documentation**: Location of the fortran binary database file defining the mesh - - -.. admonition:: Example of databases section - - .. code-block:: yaml - - databases: - mesh-database: /path/to/mesh_database.bin diff --git a/docs/parameter_documentation/header.rst b/docs/parameter_documentation/header.rst deleted file mode 100644 index 017f8a45a..000000000 --- a/docs/parameter_documentation/header.rst +++ /dev/null @@ -1,46 +0,0 @@ -Header -###### - -The header section is used for naming the run. - -Parameter definitions -===================== - -**Parameter name** : ``header`` -------------------------------- - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Define header section - -**Parameter name** : ``header.title`` -***************************************** - -**default value**: None - -**possible values**: [string] - -**documentation**: Brief name for this simulation - -**Parameter name** : ``header.description`` -********************************************* - -**default value**: None - -**possible values**: [string] - -**documentation**: Detailed description for this run. - -.. admonition:: Example header - - .. code-block:: yaml - - header: - title: Heterogeneous acoustic-elastic medium with 1 acoustic-elastic interface # name for your simulation - description: | - Material systems : Elastic domain (1), Acoustic domain (1) - Interfaces : Acoustic-elastic interface (1) - Sources : Force source (1) - Boundary conditions : Neumann BCs on all edges diff --git a/docs/parameter_documentation/index.rst b/docs/parameter_documentation/index.rst index ecac2253d..432514193 100644 --- a/docs/parameter_documentation/index.rst +++ b/docs/parameter_documentation/index.rst @@ -3,14 +3,885 @@ SPECFEM++ Parameter Documentation ================================= -The run-time behaviour of the code is controlled by YAML parameter files. The parameter file is split into a heirarchical structure of YAML nodes. -.. toctree:: - :maxdepth: 1 +On this page, we first show an example of a parameter file for SPECFEM++ and then +provide detailed documentation for each parameter in the file in a collapsible +format since a lot of the parameters are optional and only required if parent +parameters are defined. - header - simulation_setup - sources - receivers - run_setup - databases +Example parameter file +---------------------- + + +.. dropdown:: ``specfem_config.yaml`` + :open: + + .. code-block:: yaml + + parameters: + + header: + ## Header information is used for logging. It is good practice to give your simulations explicit names + title: Isotropic Elastic simulation # name for your simulation + # A detailed description for your simulation + description: | + Material systems : Elastic domain (1) + Interfaces : None + Sources : Force source (1) + Boundary conditions : Neumann BCs on all edges + + simulation-setup: + ## quadrature setup + quadrature: + quadrature-type: GLL4 + + ## Solver setup + solver: + time-marching: + time-scheme: + type: Newmark + dt: 1.1e-3 + nstep: 1600 + + simulation-mode: + forward: + writer: + seismogram: + format: "ascii" + directory: path/to/output/folder + + receivers: + stations: path/to/stations_file + angle: 0.0 + seismogram-type: + - velocity + nstep_between_samples: 1 + + ## Runtime setup + run-setup: + number-of-processors: 1 + number-of-runs: 1 + + ## databases + databases: + mesh-database: /path/to/mesh_database.bin + + ## sources + sources: path/to/sources.yaml + + + +Parameter definitions +--------------------- + + +.. dropdown:: ``parameters`` + :open: + + .. dropdown:: ``header`` + + Define header section for your simulation. This section is used for naming + the run. but has no impact on the simulation itself. + + .. dropdown:: ``title`` + + Brief name for this simulation + + :default value: None + + :possible values: [string] + + .. code-block:: yaml + :caption: Example + + title: Isotropic Elastic simulation + + .. dropdown:: ``description`` + + Detailed description for this run. + + :default value: None + + :possible values: [string] + + .. code-block:: yaml + :caption: Example + + description: | + Material systems : Elastic domain (1) + Interfaces : None + Sources : Force source (1) + Boundary conditions : Neumann BCs on all edges + + + .. dropdown:: ``simulation-setup`` + + Section to define the simulation parameters + + .. dropdown:: ``quadrature`` [optional] + + Type of quadrature used for the simulation. There are 2 ways to + define the 4th order GLL quadrature + + 1. Using predefined quadrature type + + .. code-block:: yaml + + quadrature: + quadrature-type: GLL4 + + 2. Using individual parameters + + .. code-block:: yaml + + quadrature: + alpha: 0.0 + beta: 0.0 + ngllx: 5 + ngllz: 5 + + .. dropdown:: ``quadrature-type`` [optional] + + Predefined quadrature types. + + 1. ``GLL4`` defines 4th order GLL quadrature with 5 GLL points. + 2. ``GLL7`` defines 7th order GLL quadrature with 8 GLL points. + + :default value: GLL4 + + :possible values: [GLL4, GLL7] + + + .. dropdown:: ``alpha`` + + Alpha value of the Gauss-Jacobi quadrature. For GLL quadrature alpha + = 0.0 + + :default value: None + + :possible values: [float, double] + + .. code-block:: yaml + :caption: Example + + quadrature: + alpha: 0.0 + + + .. dropdown:: ``beta`` + + Beta value of the Gauss-Jacobi quadrature. For GLL quadrature beta = + 0.0, and for GLJ quadrature beta = 1.0 + + :default value: None + + :possible values: [float, double] + + .. code-block:: yaml + :caption: Example + + quadrature: + beta: 0.0 + + + .. dropdown:: ``ngllx`` + + Number of GLL points in ``X`` dimension. + + :default value: None + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + quadrature: + ngllx: 5 + + + .. dropdown:: ``ngllz`` + + Number of GLL points in ``X`` dimension. + + :default value: None + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + quadrature: + ngllz: 5 + + + .. dropdown:: ``solver`` + + Section to define the type of solver to use for the simulation. + + .. code-block:: yaml + :caption: Example solver section + + solver: + time-marching: + time-scheme: + type: Newmark + dt: 0.001 + nstep: 1000 + t0: 0.0 + + + .. dropdown:: time-marching + + Select either a time-marching or an explicit solver. Only + time-marching solver is implemented currently. + + .. dropdown:: ``time-scheme`` + + Section to define the time scheme for the solver. + + .. dropdown:: ``type`` + + Select time scheme for the solver + + :default value: None + + :possible values: [Newmark] + + .. code-block:: yaml + :caption: Example + + time-scheme: + type: Newmark + + + .. dropdown:: ``dt`` + + Value of time step in seconds + + :default value: None + + :possible values: [float, double] + + .. code-block:: yaml + :caption: Example + + time-scheme: + dt: 0.001 + + + .. dropdown:: ``nstep`` + + Total number of time steps in the simulation. + + :default value: None + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + time-scheme: + nstep: 1000 + + + .. dropdown:: ``t0`` [optional] + + Start time of the simulation. + + :default value: 0.0 + + :possible values: [float, double] + + .. code-block:: yaml + :caption: Example + + time-scheme: + t0: 0.0 + + .. dropdown:: ``simulation-mode`` + + Defines the type of simulation to run (e.g. forward, adjoint, combined, + etc.) + + .. code-block:: yaml + :caption: Example + + simulation-mode: + forward: + ... + # or + combined: + ... + + .. note:: + + Exactly one of forward or combined simulation nodes should be + defined. + + .. dropdown:: ``forward`` + + Section to define the forward solver simulation parameters. + + .. code-block:: yaml + :caption: Example forward simulation node + + forward: + writer: + seismogram: + format: ASCII + directory: /path/to/output/folder + + wavefield: + format: HDF5 + directory: /path/to/output/folder + + display: + format: PNG + directory: /path/to/output/folder + field: displacement + simulation-field: forward + time-interval: 10 + + .. note:: + + At least one writer node should be defined in the forward simulation node. + + + .. dropdown:: ``writer`` + + Defines the outputs to be stored to disk during the forward + simulation. + + .. dropdown:: ``seismogram`` + + Seismogram writer parameters. + + .. code-block:: yaml + + writer: + seismogram: + format: ASCII + directory: /path/to/output/folder + + .. dropdown:: ``format`` [optional] + + Output format of the seismogram. + + :default value: ASCII + + :possible values: [ASCII] + + + .. dropdown:: ``directory`` [optional] + + Output folder for the seismogram. + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``wavefield`` + + Forward wavefield writer parameters. + + .. code-block:: yaml + :caption: Example + + writer: + wavefield: + format: HDF5 + directory: /path/to/output/folder + + + .. dropdown:: ``format`` [optional] + + Output format of the wavefield. + + :default value: ASCII + + :possible values: [ASCII, HDF5] + + + .. dropdown:: ``directory`` [optional] + + Output folder for the wavefield. + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``display`` + + Plot the wavefield during the forward simulation. + + .. code-block:: yaml + :caption: Example + + writer: + display: + format: PNG + directory: /path/to/output/folder + field: displacement + simulation-field: forward + time-interval: 10 + + .. dropdown:: ``format`` [optional] + + Output format for resulting plots. + + :default value: PNG + + :possible values: [PNG, JPG, on_screen] + + + .. dropdown:: ``directory`` [optional] + + Output folder for the plots (not applicable for + on_screen). + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``field`` + + Component of the wavefield to be plotted. + + :default value: None + + :possible values: [displacement, velocity, acceleration, pressure] + + + .. dropdown:: ``simulation-field`` + + Type of wavefield to be plotted. + + :default value: None + + :possible values: [forward] + + + .. dropdown:: ``time-interval`` + + Time step interval for plotting the wavefield. + + :default value: None + + :possible values: [int] + + + .. dropdown:: ``combined`` [optional] + + Combined (forward + adjoint) simulation parameters. + + .. code-block:: yaml + :caption: Example combined simulation node + + simulation-mode: + combined: + reader: + wavefield: + format: HDF5 + directory: /path/to/input/folder + + ## This example avoids writing seismograms + writer: + kernels: + format: HDF5 + directory: /path/to/output/folder + + display: + format: PNG + directory: /path/to/output/folder + field: displacement + simulation-field: adjoint + time-interval: 10 + + .. note:: + + Exactly one of forward or combined simulation nodes should + be defined. + + + .. dropdown:: ``reader`` [optional] + + Defines the inputs to be read from disk during the combined + simulation. + + .. dropdown:: ``wavefield`` + + Wavefield reader parameters. + + :default value: None + + :possible values: [YAML Node] + + + .. dropdown:: ``format`` [optional] + + Format of the wavefield to be read. + + :default value: ASCII + + :possible values: [ASCII, HDF5] + + + .. dropdown:: ``directory`` [optional] + + Folder containing the wavefield to be read. + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``writer`` [optional] + + Defines the outputs to be stored to disk during the combined + simulation. + + .. dropdown:: ``seismogram`` [optional] + + Seismogram writer parameters. + + .. dropdown:: ``format`` [optional] + + Output format of the seismogram. + + :default value: ASCII + + :possible values: [ASCII] + + + .. dropdown:: ``directory`` [optional] + + Output folder for the seismogram. + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``kernels`` + + Kernel writer parameters. + + .. dropdown:: ``format`` [optional] + + Output format of the kernels. + + :default value: ASCII + + :possible values: [ASCII, HDF5] + + + .. dropdown:: ``directory`` [optional] + + Output folder for the kernels. + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``display`` [optional] + + Plot the wavefield during the combined simulation. + + .. dropdown:: ``format`` [optional] + + Output format for resulting plots. + + :default value: PNG + + :possible values: [PNG, JPG, on_screen] + + + .. dropdown:: ``directory`` [optional] + + Output folder for the plots (not applicable for + on_screen). + + :default value: Current working directory + + :possible values: [string] + + + .. dropdown:: ``field`` + + Component of the wavefield to be plotted. + + :default value: None + + :possible values: [displacement, velocity, acceleration, pressure] + + + .. dropdown:: ``simulation-field`` + + Type of wavefield to be plotted. + + :default value: None + + :possible values: [adjoint, backward] + + + .. dropdown:: ``time-interval`` + + Time step interval for plotting the wavefield. + + :default value: None + + :possible values: [int] + + + .. dropdown:: ``receivers`` + + Parameter file section that contains the receiver information required to + calculate seismograms. + + .. code-block:: yaml + :caption: Example receivers section + + receivers: + stations: /path/to/stations_file + angle: 0.0 + seismogram-type: + - velocity + - displacement + nstep_between_samples: 1 + + .. note:: + + Please note that the ``stations_file`` is generated using SPECFEM2D mesh + generator i.e. xmeshfem2d + + .. dropdown:: ``stations`` + + Path to ``stations_file``. + + :default value: None + + :possible values: [string] + + .. code-block:: yaml + :caption: Example + + stations: /path/to/stations_file + + + .. dropdown:: ``angle`` + + Angle to rotate components at receivers + + :default value: None + + :possible values: [float] + + .. code-block:: yaml + :caption: Example + + angle: 0.0 + + + .. dropdown:: ``seismogram-type`` + + Type of seismograms to be written. + + :default value: None + + :possible values: [YAML list] + + .. code-block:: yaml + :caption: Example + + seismogram-type: + - velocity + - displacement + + .. rst-class:: center-table + + +-------------------+---------------------------------------+-------------------------------------+ + | Seismogram | SPECFEM Par_file ``seismotype`` value | ``receivers.seismogram-type`` value | + +===================+=======================================+=====================================+ + | Displacement | 1 | ``displacement`` | + +-------------------+---------------------------------------+-------------------------------------+ + | Velocity | 2 | ``velocity`` | + +-------------------+---------------------------------------+-------------------------------------+ + | Acceleration | 3 | ``acceleration`` | + +-------------------+---------------------------------------+-------------------------------------+ + | Pressure | 4 | ``pressure`` | + +-------------------+---------------------------------------+-------------------------------------+ + | Displacement Curl | 5 | ✘ Unsupported | + +-------------------+---------------------------------------+-------------------------------------+ + | Fluid Potential | 6 | ✘ Unsupported | + +-------------------+---------------------------------------+-------------------------------------+ + + + + .. dropdown:: ``nstep_between_samples`` + + Number of time steps between sampling the wavefield at station locations + for writing seismogram. + + :default value: None + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + nstep_between_samples: 1 + + + + .. dropdown:: ``run-setup`` + + Define run-time configuration for your simulation. + + .. code-block:: yaml + :caption: Example run-setup section + + run-setup: + number-of-processors: 1 + number-of-runs: 1 + + .. dropdown:: ``number-of-processors`` + + Number of MPI processes used in the simulation. MPI version is not + enabled in this version of the package. number-of-processors == 1 + + :default value: 1 + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + number-of-processors: 1 + + + .. dropdown:: ``number-of-runs`` + + Number of runs in this simulation. Only single run implemented in this + version of the package. number-of-runs == 1 + + :default value: 1 + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + number-of-runs: 1 + + + + .. dropdown:: ``databases`` + + The databases section defines the location of files to be read by the + solver. + + .. code-block:: yaml + :caption: Example of databases section + + databases: + mesh-database: /path/to/mesh_database.bin + + + .. _database-file-parameter: + + .. dropdown:: ``mesh-database`` + :open: + + Location of the fortran binary database file defining the mesh + + :default value: None + + :possible values: [string] + + .. code-block:: yaml + :caption: Example + + mesh-database: /path/to/mesh_database.bin + + + .. dropdown:: ``sources`` + + Define sources + + :default value: None + + :possible values: [string, YAML Node] + + .. admonition:: Example sources section + + The sources is a path to a YAML file. + + .. code-block:: yaml + + sources: path/to/sources.yaml + + The sources section is a YAML node that contains the source information + + .. code-block:: yaml + + sources: + number-of-sources: 1 + sources: + - force: + x : 2500.0 + z : 2500.0 + source_surf: false + angle : 0.0 + vx : 0.0 + vz : 0.0 + Ricker: + factor: 1e10 + tshift: 0.0 + f0: 10.0 + + .. note:: + + The parameters below are only relevant if the sources section is + defined as a YAML node. + + .. dropdown:: ``number-of-sources`` + + Number of sources in the simulation + + :default value: None + + :possible values: [int] + + .. code-block:: yaml + :caption: Example + + number-of-sources: 1 + + + .. dropdown:: ``sources`` + + List of sources + + :default value: None + + :possible values: [YAML Node] + + .. code-block:: yaml + :caption: Example + + sources: + - force: + x : 2500.0 + z : 2500.0 + source_surf: false + angle : 0.0 + vx : 0.0 + vz : 0.0 + Ricker: + factor: 1e10 + tshift: 0.0 + f0: 10.0 diff --git a/docs/parameter_documentation/receivers.rst b/docs/parameter_documentation/receivers.rst deleted file mode 100644 index fa2256a21..000000000 --- a/docs/parameter_documentation/receivers.rst +++ /dev/null @@ -1,90 +0,0 @@ -receivers -########## - -Receivers section defines receiver information required to calculate seismograms. - -.. note:: - - Please note that the ``stations_file`` is generated using SPECFEM2D mesh generator i.e. xmeshfem2d - -**Parameter Name** : ``receivers`` ------------------------------------ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : receiver information required to calculate seismograms. - -**Parameter Name** : ``receivers.stationse`` -****************************************************** - -**default value** : None - -**possible values** : [string] - -**documentation** : Path to ``stations_file`` - -**Parameter Name** : ``receivers.angle`` -****************************************************** - -**default value** : None - -**possible values** : [float] - -**documentation** : Angle to rotate components at receivers - -**Parameter Name** : ``receivers.seismogram-type`` -************************************************** - -**default value** : None - -**possible values** : [YAML list] - -**documentation** : Type of seismograms to be written. - - .. rst-class:: center-table - - +-------------------+---------------------------------------+-------------------------------------+ - | Seismogram | SPECFEM Par_file ``seismotype`` value | ``receivers.seismogram-type`` value | - +===================+=======================================+=====================================+ - | Displacement | 1 | ``displacement`` | - +-------------------+---------------------------------------+-------------------------------------+ - | Velocity | 2 | ``velocity`` | - +-------------------+---------------------------------------+-------------------------------------+ - | Acceleration | 3 | ``acceleration`` | - +-------------------+---------------------------------------+-------------------------------------+ - | Pressure | 4 | ``pressure`` | - +-------------------+---------------------------------------+-------------------------------------+ - | Displacement Curl | 5 | ✘ Unsupported | - +-------------------+---------------------------------------+-------------------------------------+ - | Fluid Potential | 6 | ✘ Unsupported | - +-------------------+---------------------------------------+-------------------------------------+ - - -.. code-block:: yaml - - seismogram-type: - - velocity - - displacement - -**Parameter Name** : ``receivers.nstep_between_samples`` -******************************************************** - -**default value** : None - -**possible values** : [int] - -**documentation** : Number of time steps between sampling the wavefield at station locations for writing seismogram. - -.. admonition:: Example receivers section - - .. code-block:: yaml - - receivers: - stations: /path/to/stations_file - angle: 0.0 - seismogram-type: - - velocity - - displacement - nstep_between_samples: 1 diff --git a/docs/parameter_documentation/run_setup.rst b/docs/parameter_documentation/run_setup.rst deleted file mode 100644 index 30a14b9a6..000000000 --- a/docs/parameter_documentation/run_setup.rst +++ /dev/null @@ -1,42 +0,0 @@ -Runtime Setup -############# - -Runtime setup section defines the run-time setup of the simulation. If this section is not defined, the simulation will be a serial simulation with a single run. - -Parameter definitions -===================== - -**Parameter Name** : ``run-setup`` [optional] ------------------------------------ - -**default value** : number-of-processors: 1, number-of-runs: 1 - -**possible values** : [YAML Node] - -**documentation** : Define run-time configuration for your simulation - -**Parameter Name** : ``run-setup.number-of-processors`` [optional] -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -**default value** : 1 - -**possible values** : [int] - -**documentation** : Number of MPI processes used in the simulation. MPI version is not enabled in this version of the package. number-of-processors == 1 - -**Parameter Name** : ``run-setup.number-of-runs`` [optional] -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -**default value** : 1 - -**possible values** : [int] - -**documentation** : Number of runs in this simulation. Only single run implemented in this version of the package. number-of-runs == 1 - -.. admonition:: Example run-setup section - - .. code-block:: yaml - - run-setup: - number-of-processors: 1 - number-of-runs: 1 diff --git a/docs/parameter_documentation/seismogram_setup.rst b/docs/parameter_documentation/seismogram_setup.rst deleted file mode 100644 index 3c6771ae3..000000000 --- a/docs/parameter_documentation/seismogram_setup.rst +++ /dev/null @@ -1,37 +0,0 @@ -Seismograms -########### - -Seismograms section defines the output seismograms format at predefined stations. When not-defined seismograms will be calculated but not written to a file. - - -Parameter definitions -======================= - -**Parameter Name** : ``seismogram`` ------------------------------------- - -**default value** : NULL - -**possible values** : [YAML Node] - -**documentation** : Define seismogram output configuration - -**Parameter Name** : ``seismogram.seismogram-format`` -****************************************************** - -**default value** : None - -**possible values** : [ascii] - -**documentation** : Type of seismogram format to be written. - -1. ascii - :ref:`ASCII` writes calculated seismogram values to seismogram files in string format. - -**Parameter Name** : ``seismogram.output-folder`` -****************************************************** - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Path to output folder where the seismograms will be saved. diff --git a/docs/parameter_documentation/simulation_setup.rst b/docs/parameter_documentation/simulation_setup.rst deleted file mode 100644 index 4e8ca46ee..000000000 --- a/docs/parameter_documentation/simulation_setup.rst +++ /dev/null @@ -1,513 +0,0 @@ -Simulation Setup -################ - -Simulation setup defines the run-time behaviour of the simulation - -Parameter definitions -===================== - -**Parameter Name** : ``simulation-setup`` ------------------------------------------ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Simulation setup parameters - -**Parameter Name** : ``simulation-setup.quadrature`` [optional] ----------------------------------------------- - -**default value** : 4th order GLL quadrature with 5 GLL points - -**possible values**: [YAML Node] - -**documentation** : Type of quadrature used for the simulation. - -**Parameter Name** : ``simulation-setup.quadrature.alpha`` [optional] -**************************************************** - -**default value** : None - -**possible values** : [float, double] - -**documentation** : Alpha value of the Gauss-Jacobi quadrature. For GLL quadrature alpha = 0.0 - -**Parameter Name** : ``simulation-setup.quadrature.beta`` [optional] -*************************************************** - -**default value** : None - -**possible values** : [float, double] - -**documentation** : Beta value of the Gauss-Jacobi quadrature. For GLL quadrature beta = 0.0, and for GLJ quadrature beta = 1.0 - -**Parameter Name** : ``simulation-setup.quadrature.ngllx`` [optional] -**************************************************** - -**default value** : None - -**possible values** : [int] - -**documentation** : Number of GLL points in X-dimension - -**Parameter Name** : ``simulation-setup.quadrature.ngllz`` [optional] -***************************************************** - -**default value** : None - -**possible values** : [int] - -**documentation** : Number of GLL points in Z-dimension - -**Parameter Name** : ``simulation-setup.quadrature.quadrature-type`` [optional] -************************************************************** - -**default value** : GLL4 - -**possible values** : [GLL4, GLL7] - -**decumentation** : Predefined quadrature types. - -1. ``GLL4`` defines 4th order GLL quadrature with 5 GLL points. -2. ``GLL7`` defines 7th order GLL quadrature with 8 GLL points. - -.. admonition:: Example for defining 4th order GLL quadrature - - There are 2 ways to define the 4th order GLL quadrature - - 1. Using predefined quadrature type - - .. code-block:: yaml - - quadrature: - quadrature-type: GLL4 - - 2. Using individual parameters - - .. code-block:: yaml - - quadrature: - alpha: 0.0 - beta: 0.0 - ngllx: 5 - ngllz: 5 - -**Parameter Name** : ``simulation-setup.solver`` -------------------------------- - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Type of solver to use for the simulation. - -**Parameter Name** : ``simulation-setup.solver.time-marching`` -********************************************** - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Select either a time-marching or an explicit solver. Only time-marching solver is implemented currently. - - -**Parameter Name** : ``simulation-setup.solver.time-marching.time-scheme.type`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [Newmark] - -**documentation** : Select time scheme for the solver - -**Parameter Name** : ``simulation-setup.solver.time-marching.time-scheme.dt`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [float, double] - -**documentation** : Value of time step in seconds - -**Parameter Name** : ``simulation-setup.solver.time-marching.time-scheme.nstep`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [int] - -**documentation** : Total number of time steps in the simulation - -**Parameter Name** : ``simulation-setup.solver.time-marching.time-scheme.t0`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : 0.0 - -**possible values** : [float, double] - -**documentation** : Start time of the simulation - -.. admonition:: Example for defining time-marching Newmark solver - - .. code-block:: yaml - - solver: - time-marching: - time-scheme: - type: Newmark - dt: 0.001 - nstep: 1000 - t0: 0.0 - -**Parameter Name** : ``simulation-setup.simulation-mode`` ---------------------------------------------------------- - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Defines the type of simulation to run (e.g. forward, adjoint, combined, etc.) - -**Parameter Name** : ``simulation-setup.simulation-mode.forward`` [optional] -***************************************************************************** - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Forward simulation parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Defines the outputs to be stored to disk during the forward simulation - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.seismogram`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Seismogram writer parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.seismogram.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : ASCII - -**possible values** : [ASCII] - -**documentation** : Output format of the seismogram - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.seismogram.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Output folder for the seismogram - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.wavefield`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Forward wavefield writer parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.wavefield.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : ASCII - -**possible values** : [ASCII, HDF5] - -**documentation** : Output format of the wavefield - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.wavefield.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Output folder for the wavefield - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display`` [optional] -******************************************************************************************* - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Plot the wavefield during the forward simulation - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : PNG - -**possible values** : [PNG, JPG, on_screen] - -**documentation** : Output format for resulting plots - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Output folder for the plots (not applicable for on_screen) - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.field`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [displacement, velocity, acceleration, pressure] - -**documentation** : Component of the wavefield to be plotted - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.simulation-field`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [forward] - -**documentation** : Type of wavefield to be plotted - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.time-interval`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [int] - -**documentation** : Time step interval for plotting the wavefield - -.. admonition:: Example for defining a forward simulation node - - .. code-block:: yaml - - simulation-mode: - forward: - writer: - seismogram: - format: ASCII - directory: /path/to/output/folder - - wavefield: - format: HDF5 - directory: /path/to/output/folder - - display: - format: PNG - directory: /path/to/output/folder - field: displacement - simulation-field: forward - time-interval: 10 - - -.. Note:: - - Atlease one writer node should be defined in the forward simulation node. - -**Parameter Name** : ``simulation-setup.simulation-mode.combined`` [optional] -***************************************************************************** - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Combined (forward + adjoint) simulation parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.reader`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Defines the inputs to be read from disk during the combined simulation - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.reader.wavefield`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Wavefield reader parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.reader.wavefield.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : ASCII - -**possible values** : [ASCII, HDF5] - -**documentation** : Format of the wavefield to be read - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.reader.wavefield.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Folder containing the wavefield to be read - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Defines the outputs to be stored to disk during the combined simulation - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer.seismogram`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Seismogram writer parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer.seismogram.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : ASCII - -**possible values** : [ASCII] - -**documentation** : Output format of the seismogram - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer.seismogram.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Output folder for the seismogram - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer.kernels`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Kernel writer parameters - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer.kernels.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : ASCII - -**possible values** : [ASCII, HDF5] - -**documentation** : Output format of the kernels - -**Parameter Name** : ``simulation-setup.simulation-mode.combined.writer.kernels.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Output folder for the kernels - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [YAML Node] - -**documentation** : Plot the wavefield during the forward simulation - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.format`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : PNG - -**possible values** : [PNG, JPG, on_screen] - -**documentation** : Output format for resulting plots - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.directory`` [optional] -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : Current working directory - -**possible values** : [string] - -**documentation** : Output folder for the plots (not applicable for on_screen) - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.field`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [displacement, velocity, acceleration, pressure] - -**documentation** : Component of the wavefield to be plotted - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.simulation-field`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [adjoint, backward] - -**documentation** : Type of wavefield to be plotted - -**Parameter Name** : ``simulation-setup.simulation-mode.forward.writer.display.time-interval`` -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -**default value** : None - -**possible values** : [int] - -**documentation** : Time step interval for plotting the wavefield - -.. admonition:: Example for defining a combined simulation node - - .. code-block:: yaml - - simulation-mode: - combined: - reader: - wavefield: - format: HDF5 - directory: /path/to/input/folder - - ## This example avoids writing seismograms - writer: - kernels: - format: HDF5 - directory: /path/to/output/folder - - display: - format: PNG - directory: /path/to/output/folder - field: displacement - simulation-field: adjoint - time-interval: 10 - -.. Note:: - - Exactly one of forward or combined simulation nodes should be defined. diff --git a/docs/parameter_documentation/sources.rst b/docs/parameter_documentation/sources.rst deleted file mode 100644 index a74fd9eee..000000000 --- a/docs/parameter_documentation/sources.rst +++ /dev/null @@ -1,15 +0,0 @@ -**Parameter name** : ``sources`` -****************************************************** - -**default value**: None - -**possible values**: [string] - -**documentation**: Location of source file (yaml) defining the location of sources - -.. admonition:: Example of databases section - - .. code-block:: yaml - - sources: - sources: path/to/sources.yaml diff --git a/docs/requirements.txt b/docs/requirements.txt index b55aacb5b..a63d16f77 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,4 +1,6 @@ breathe furo +sphinx==8.0.0 sphinx-copybutton sphinx-sitemap +sphinx_design diff --git a/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 b/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 index 82435a529..9a334ad8f 100644 --- a/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 +++ b/fortran/meshfem3d/generate_databases/save_arrays_solver.F90 @@ -75,13 +75,94 @@ subroutine save_global_arrays_with_components(nspec, array) end subroutine save_global_arrays_with_components - end module save_arrays_module + subroutine save_boundary_arrays(num_faces, spec_boundary, ijk_boundary, & + jacobian2Dw_boundary, normal_boundary) + + use constants, only: IOUT, CUSTOM_REAL + implicit none + + integer, intent(in) :: num_faces + integer :: iface + integer, dimension(:), intent(in) :: spec_boundary + integer, dimension(:,:,:), intent(in) :: ijk_boundary + real(kind=CUSTOM_REAL), dimension(:,:), intent(in) :: jacobian2Dw_boundary + real(kind=CUSTOM_REAL), dimension(:,:,:), intent(in) :: normal_boundary + + ! Save the spec_boundary array + write(iout) spec_boundary + + ! save the ijk_boundary array + do iface = 1, num_faces + write(iout) ijk_boundary(:, :, iface) + end do + + ! save the jacobian2dw_boundary array + do iface = 1, num_faces + write(iout) jacobian2dw_boundary(:, iface) + end do + + ! save the normal_boundary array + do iface = 1, num_faces + write(iout) normal_boundary(:, :, iface) + end do + + 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_ispec_is_arrays(nspec, ispec_is_acoustic, ispec_is_elastic, ispec_is_poroelastic) + + use constants, only: IOUT + + 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 + 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 + 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_ispec_is_arrays + +end module save_arrays_module ! for external mesh subroutine save_arrays_solver_mesh() use constants, only: IMAIN,IOUT,myrank - use save_arrays_module, only: save_global_arrays, save_global_arrays_with_components + use save_arrays_module, only: & + save_global_arrays, & + save_global_arrays_with_components, & + save_ispec_is_arrays, & + save_boundary_arrays use shared_parameters, only: ACOUSTIC_SIMULATION, ELASTIC_SIMULATION, POROELASTIC_SIMULATION, & APPROXIMATE_OCEAN_LOAD, SAVE_MESH_FILES, ANISOTROPY @@ -205,25 +286,26 @@ subroutine save_arrays_solver_mesh() call save_global_arrays(nspec, gammazstore) call save_global_arrays(nspec, jacobianstore) - ! write(IOUT) xixstore - ! write(IOUT) xiystore - ! write(IOUT) xizstore - ! write(IOUT) etaxstore - ! write(IOUT) etaystore - ! write(IOUT) etazstore - ! write(IOUT) gammaxstore - ! write(IOUT) gammaystore - ! write(IOUT) gammazstore - ! write(IOUT) jacobianstore + ! write test value + itest = 10000 + write(IOUT) itest call save_global_arrays(nspec, kappastore) call save_global_arrays(nspec, mustore) - ! write(IOUT) kappastore - ! write(IOUT) mustore - write(IOUT) - 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 @@ -236,10 +318,11 @@ subroutine save_arrays_solver_mesh() ! this array is needed for acoustic simulations but also for elastic simulations with CPML, ! thus we allocate it and read it in all cases (whether the simulation is acoustic, elastic, or acoustic/elastic) - call save_global_arrays(nspec, rhostore) - ! write(IOUT) rhostore + ! write test value + itest = 9998 + write(IOUT) itest ! elastic if (ELASTIC_SIMULATION) then @@ -250,10 +333,12 @@ subroutine save_arrays_solver_mesh() ! Stacey call save_global_arrays(nspec, rho_vp) call save_global_arrays(nspec, rho_vs) - ! write(IOUT) rho_vp - ! write(IOUT) rho_vs endif + ! Write a test value + itest = 9997 + write(IOUT) itest + ! poroelastic if (POROELASTIC_SIMULATION) then write(IOUT) rmass_solid_poroelastic @@ -269,17 +354,12 @@ subroutine save_arrays_solver_mesh() call save_global_arrays(nspec, rho_vpII) call save_global_arrays(nspec, rho_vsI) - ! write(IOUT) rhoarraystore - ! write(IOUT) kappaarraystore - ! write(IOUT) etastore - ! write(IOUT) tortstore - ! write(IOUT) permstore - ! write(IOUT) phistore - ! write(IOUT) rho_vpI - ! write(IOUT) rho_vpII - ! write(IOUT) rho_vsI endif + ! write test value + itest = 9996 + write(IOUT) itest + ! @Lucas & @Congyue need to uncomment this when implementing PML ! C-PML absorbing boundary conditions @@ -318,10 +398,14 @@ subroutine save_arrays_solver_mesh() ! absorbing boundary surface write(IOUT) num_abs_boundary_faces if (num_abs_boundary_faces > 0) then - write(IOUT) abs_boundary_ispec - write(IOUT) abs_boundary_ijk - write(IOUT) abs_boundary_jacobian2Dw - write(IOUT) abs_boundary_normal + + ! Saves the arrays with num faces in the first dimension + call save_boundary_arrays(num_abs_boundary_faces, & + abs_boundary_ispec, & + abs_boundary_ijk, & + abs_boundary_jacobian2Dw, & + abs_boundary_normal) + if (STACEY_ABSORBING_CONDITIONS .and. (.not. PML_CONDITIONS)) then ! store mass matrix contributions if (ELASTIC_SIMULATION) then @@ -336,7 +420,7 @@ subroutine save_arrays_solver_mesh() endif ! stamp for checking i/o so far - itest = 9998 + itest = 9995 write(IOUT) itest ! boundaries @@ -357,39 +441,58 @@ subroutine save_arrays_solver_mesh() ! free surface write(IOUT) num_free_surface_faces if (num_free_surface_faces > 0) then - write(IOUT) free_surface_ispec - write(IOUT) free_surface_ijk - write(IOUT) free_surface_jacobian2Dw - write(IOUT) free_surface_normal + + ! Saves the arrays with num faces in the first dimension + call save_boundary_arrays(num_free_surface_faces, & + free_surface_ispec, & + free_surface_ijk, & + free_surface_jacobian2Dw, & + free_surface_normal) endif ! acoustic-elastic coupling surface write(IOUT) num_coupling_ac_el_faces if (num_coupling_ac_el_faces > 0) then - write(IOUT) coupling_ac_el_ispec - write(IOUT) coupling_ac_el_ijk - write(IOUT) coupling_ac_el_jacobian2Dw - write(IOUT) coupling_ac_el_normal + + ! Saves the arrays with num faces in the first dimension + call save_boundary_arrays(num_coupling_ac_el_faces, & + coupling_ac_el_ispec, & + coupling_ac_el_ijk, & + coupling_ac_el_jacobian2Dw, & + coupling_ac_el_normal) + endif ! acoustic-poroelastic coupling surface write(IOUT) num_coupling_ac_po_faces if (num_coupling_ac_po_faces > 0) then - write(IOUT) coupling_ac_po_ispec - write(IOUT) coupling_ac_po_ijk - write(IOUT) coupling_ac_po_jacobian2Dw - write(IOUT) coupling_ac_po_normal + + ! Saves the arrays with num faces in the first dimension + call save_boundary_arrays(num_coupling_ac_po_faces, & + coupling_ac_po_ispec, & + coupling_ac_po_ijk, & + coupling_ac_po_jacobian2Dw, & + coupling_ac_po_normal) endif ! elastic-poroelastic coupling surface write(IOUT) num_coupling_el_po_faces if (num_coupling_el_po_faces > 0) then - write(IOUT) coupling_el_po_ispec - write(IOUT) coupling_po_el_ispec - write(IOUT) coupling_el_po_ijk - write(IOUT) coupling_po_el_ijk - write(IOUT) coupling_el_po_jacobian2Dw - write(IOUT) coupling_el_po_normal + + ! Saves the arrays with num faces in the first dimension + call save_boundary_arrays(num_coupling_el_po_faces, & + coupling_el_po_ispec, & + coupling_el_po_ijk, & + coupling_el_po_jacobian2Dw, & + coupling_el_po_normal) + + ! Saves the arrays with num faces in the first dimension + call save_boundary_arrays(num_coupling_el_po_faces, & + coupling_po_el_ispec, & + coupling_po_el_ijk, & + coupling_el_po_jacobian2Dw, & + coupling_el_po_normal) + endif ! stamp for checking i/o diff --git a/fortran/meshfem3d/generate_databases/save_parameters.f90 b/fortran/meshfem3d/generate_databases/save_parameters.f90 index 5c8397a40..5358f2f4f 100644 --- a/fortran/meshfem3d/generate_databases/save_parameters.f90 +++ b/fortran/meshfem3d/generate_databases/save_parameters.f90 @@ -8,7 +8,7 @@ subroutine save_parameters() IOUT, myrank, sizeprocs, & NDIM, NSPEC_AB, NGLOB_AB, & NGLLX, NGLLY, NGLLZ, nfaces_surface, num_neighbors_all, NGLLSQUARE, & - nspec2D_xmin, nspec2D_xmin, nspec2D_ymin, nspec2D_ymin, & + nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, & NSPEC2D_BOTTOM, NSPEC2D_TOP, & ACOUSTIC_SIMULATION, ELASTIC_SIMULATION, POROELASTIC_SIMULATION, & ANISOTROPY, & diff --git a/include/IO/fortranio/fortran_io.tpp b/include/IO/fortranio/fortran_io.tpp index 929f13ac2..3007a824c 100644 --- a/include/IO/fortranio/fortran_io.tpp +++ b/include/IO/fortranio/fortran_io.tpp @@ -23,6 +23,14 @@ void fortran_read_value(std::vector *value, return; } + +// Template specialization for std::vector +template <> +void fortran_read_value(std::vector *value, + std::ifstream &stream, + int &buffer_length); + + template void fortran_IO(std::ifstream &stream, int &buffer_length, T *value, Args... values) { @@ -41,7 +49,16 @@ void fortran_read_line(std::ifstream &stream, Args... values) { } stream.read(reinterpret_cast(&buffer_length), fint); - specfem::IO::fortran_IO(stream, buffer_length, values...); + + try { + specfem::IO::fortran_IO(stream, buffer_length, values...); + } catch (const std::exception &e) { + std::ostringstream error_message; + error_message << "Error reading fortran line with buffer length: " + << buffer_length << "\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } stream.read(reinterpret_cast(&buffer_length), fint); return; diff --git a/include/IO/fortranio/interface.hpp b/include/IO/fortranio/interface.hpp index 3cc400dce..74c0a4e59 100644 --- a/include/IO/fortranio/interface.hpp +++ b/include/IO/fortranio/interface.hpp @@ -10,10 +10,23 @@ namespace IO { /** * @brief Read a line from fortran unformatted binary file * - * @tparam Args Argument can be of the type bool, int, type_real, string, - * vector + * This function is the core of the fortran binary file reading. It reads a line + * from the file and assigns the values to the variables passed as arguments. + * + * @tparam Args Argument can be of the type `bool`, `int`, `type_real`, + * `string`, or `vector` * @param stream An open file stream. * @param values Comma separated list of variable addresses to be read. + * @throws std::runtime_error if an error occurs while reading the line + * + * @code{.cpp} // Example of how to use this function + * int value1, value2; + * specfem::IO::fortran_read_line(stream, &value1, &value2); + * + * // Example of how to use this function with a vector + * std::vector values(10); + * specfem::IO::fortran_read_line(stream, &values); + * @endcode */ template void fortran_read_line(std::ifstream &stream, Args... values); diff --git a/include/IO/interface.hpp b/include/IO/interface.hpp index c0d6da19e..ffa46a423 100644 --- a/include/IO/interface.hpp +++ b/include/IO/interface.hpp @@ -18,16 +18,27 @@ namespace IO { * * @param filename Fortran binary database filename * @param mpi pointer to MPI object to manage communication - * @return specfem::mesh::mesh Specfem mesh object + * @return specfem::mesh::mesh Specfem mesh object for dimension type dim2 * */ specfem::mesh::mesh read_2d_mesh(const std::string filename, const specfem::MPI::MPI *mpi); +/** + * @brief Construct a 3D mesh object from a Fortran binary database file + * + * @param mesh_parameters_file Mesh parameters file + * @param mesh_databases_file Mesh databases file + * @param mpi pointer to MPI object to manage communication + * @return specfem::mesh::mesh + * Specfem mesh object for dimension type dim3 + * + */ specfem::mesh::mesh read_3d_mesh(const std::string mesh_parameters_file, const std::string mesh_databases_file, const specfem::MPI::MPI *mpi); + /** * @brief Read station file * @@ -73,8 +84,12 @@ read_receivers(const YAML::Node &stations, const type_real angle); * Parse source specification file written in yaml format and create a vector of * specfem::source::source * object * - * @param sources_file Name of the yml file - * @param mpi Pointer to specfem MPI object + * @param sources_file Name of the yaml file + * @param nsteps Number of time steps + * @param user_t0 User defined t0 + * @param dt Time step + * @param simulation_type Type of simulation + * * @return std::vector vector of instantiated source * objects */ @@ -90,7 +105,10 @@ read_sources(const std::string sources_file, const int nsteps, * specfem::source::source * object * * @param yaml YAML node containing source information - * @param mpi Pointer to specfem MPI object + * @param nsteps Number of time steps + * @param user_t0 User defined t0 + * @param dt Time step + * @param simulation_type Type of simulation * @return std::vector vector of instantiated source * objects */ diff --git a/include/IO/mesh/impl/fortran/dim2/read_material_properties.hpp b/include/IO/mesh/impl/fortran/dim2/read_material_properties.hpp index 4aed4f750..f084fa960 100644 --- a/include/IO/mesh/impl/fortran/dim2/read_material_properties.hpp +++ b/include/IO/mesh/impl/fortran/dim2/read_material_properties.hpp @@ -1,5 +1,6 @@ #pragma once +#include "enumerations/dimension.hpp" #include "mesh/mesh.hpp" #include "specfem_mpi/interface.hpp" #include @@ -25,9 +26,11 @@ namespace dim2 { * from the database file */ -specfem::mesh::materials read_material_properties( - std::ifstream &stream, const int numat, const int nspec, - const specfem::kokkos::HostView2d knods, const specfem::MPI::MPI *mpi); +specfem::mesh::materials +read_material_properties(std::ifstream &stream, const int numat, + const int nspec, + const specfem::kokkos::HostView2d knods, + const specfem::MPI::MPI *mpi); } // namespace dim2 } // namespace fortran diff --git a/include/IO/mesh/impl/fortran/dim3/interface.hpp b/include/IO/mesh/impl/fortran/dim3/interface.hpp index d965fdb7c..f9e77ba17 100644 --- a/include/IO/mesh/impl/fortran/dim3/interface.hpp +++ b/include/IO/mesh/impl/fortran/dim3/interface.hpp @@ -11,7 +11,7 @@ namespace impl { namespace fortran { namespace dim3 { -/* +/** * @brief Read paramters from 3D mesh database * * @param stream Input stream @@ -22,7 +22,7 @@ namespace dim3 { specfem::mesh::parameters read_mesh_parameters(std::ifstream &stream, const specfem::MPI::MPI *mpi); -/* +/** * @brief Read mapping from 3D mesh database * * @param stream Input stream @@ -33,7 +33,19 @@ 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 * * @param stream Input stream @@ -45,7 +57,7 @@ void read_xyz( specfem::mesh::coordinates &coordinates, const specfem::MPI::MPI *mpi); -/* +/** * @brief Read Jacobian from 3D mesh database * * @param stream Input stream @@ -58,25 +70,87 @@ void read_partial_derivatives( &partial_derivatives, const specfem::MPI::MPI *mpi); -template using View1D = Kokkos::View; - -template using View4D = Kokkos::View; - -template using View5D = Kokkos::View; - -template void read_array(std::ifstream &stream, View1D &array); - -template void read_array(std::ifstream &stream, View4D &array); +/** + * @brief Read array from 3D mesh database + * + * This function is really used to read any size array from the Fortran + * 3D mesh database. It works with 1-5D arrays, and assumes that the first + * dimension is always the last dimension of the fortran array. And unwraps + * the rest of the dimensions (which are written in Fortran order) to the + * Kokkos::View. + * + * @param stream Input stream + * @param array Array to read + * @tparam ViewType Kokkos::View type + * @throws std::runtime_error if an error occurs while reading the array + * + * @code{.cpp} + * // Example of how to use this function + * Kokkos::View array("array", 10); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, array); + * @endcode + */ +template +void read_array(std::ifstream &stream, ViewType &array); -template void read_array(std::ifstream &stream, View5D &array); +/** + * @brief Read index array from 3D mesh database, subtracts 1 from each value + * to convert from Fortran to C indexing + * + * This function is really used to read any size index array from the Fortran + * 3D mesh database. It works with 1-4D arrays, and assumes that the first + * dimension is always the last dimension of the fortran array. And unwraps + * the rest of the dimensions (which are written in Fortran order) to the + * Kokkos::View. It also subtracts 1 from each value to convert from Fortran + * to C indexing. + * + * @param stream Input stream + * @param array Index array to read + * @tparam ViewType Kokkos::View type + * @throws std::runtime_error if an error occurs while reading the array + * + * @code{.cpp} + * // Example of how to use this function + * Kokkos::View array("array", 10); + * specfem::IO::mesh::impl::fortran::dim3::read_index_array(stream, array); + * @endcode + */ +template +void read_index_array(std::ifstream &stream, ViewType &array); -// Read index array will subtract 1 from each value when reading to account for -// Fortran 1-based indexing -template -void read_index_array(std::ifstream &stream, View1D &array); +/** + * @brief Read single test value from 3D mesh database and throw error if + * value is not as expected + * + * @param stream Input stream + * @param test_value Test value to read + * @throws std::runtime_error if the test value is not as expected + * @code{.cpp} + * // Example of how to use this function + * int test_value; + * specfem::IO::mesh::impl::fortran::dim3::check_read_test_value(stream, + * test_value); + * @endcode + */ +void check_read_test_value(std::ifstream &stream, int test_value); -template -void read_index_array(std::ifstream &stream, View4D &array); +/** + * @brief Check values and throw error if they are not as expected + * + * @param message Message to print if values are not as expected + * @param value Value to check + * @param expected Expected value + * + * @throws std::runtime_error if the value is not as expected + * + * @code{.cpp} + * // Example of how to use this function + * specfem::IO::mesh::impl::fortran::dim3::check_values("message", value, + * expected); + * @endcode + * + */ +void check_values(std::string message, int value, int expected); } // namespace dim3 } // namespace fortran diff --git a/include/IO/mesh/impl/fortran/dim3/utilities.hpp b/include/IO/mesh/impl/fortran/dim3/utilities.hpp index 36dc96e88..590229565 100644 --- a/include/IO/mesh/impl/fortran/dim3/utilities.hpp +++ b/include/IO/mesh/impl/fortran/dim3/utilities.hpp @@ -1,173 +1,403 @@ +#pragma once #include "IO/fortranio/interface.hpp" #include "IO/mesh/impl/fortran/dim3/interface.hpp" #include "enumerations/dimension.hpp" #include "mesh/mesh.hpp" #include "specfem_mpi/interface.hpp" +#include "specfem_setup.hpp" #include -template using View1D = Kokkos::View; +template +void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream, + ViewType &array) { -template using View4D = Kokkos::View; + // Get the value_type and rank of the ViewType + using value_type = typename ViewType::value_type; + constexpr int rank = ViewType::rank; -template using View5D = Kokkos::View; + // 1D array implementation + if constexpr (rank == 1) { -template -void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream, - View1D &array) { - const int n = array.extent(0); - std::vector dummy_T(n, -999999); + const int n = array.extent(0); - try { - // Read into dummy vector - specfem::IO::fortran_read_line(stream, &dummy_T); - // Assign to KokkosView - for (int i = 0; i < n; i++) { - array(i) = dummy_T[i]; + std::vector dummy_T(n, -99999); + + try { + // Read into dummy vector + specfem::IO::fortran_read_line(stream, &dummy_T); + // Assign to KokkosView + for (int i = 0; i < n; i++) { + array(i) = dummy_T[i]; + } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 1D array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); } - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading array from database file:\n" - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); } -} + // 2D array implementation + else if constexpr (rank == 2) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); -template -void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream, - View4D &array) { - const int nspec = array.extent(0); - const int ngllx = array.extent(1); - const int nglly = array.extent(2); - const int ngllz = array.extent(3); + std::vector dummy_T(n1, -99999); - std::vector dummy_T(ngllx * nglly * ngllz, -999999); + try { + // Read into dummy vector + // Assign to KokkosView + for (int i = 0; i < n0; i++) { + specfem::IO::fortran_read_line(stream, &dummy_T); + int counter = 0; + for (int j = 0; j < n1; j++) { + array(i, j) = dummy_T[counter]; + counter++; + } + } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 2D array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } + } + // 3D array implementation + else if constexpr (rank == 3) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); + const int n2 = array.extent(2); - try { - for (int ispec = 0; ispec < nspec; ispec++) { - // Read into dummy vector for each ispec - specfem::IO::fortran_read_line(stream, &dummy_T); + std::vector dummy_T(n1 * n2, -99999); - // Assign to KokkosView - int counter = 0; - for (int igllx = 0; igllx < ngllx; igllx++) { - for (int iglly = 0; iglly < nglly; iglly++) { - for (int igllz = 0; igllz < ngllz; igllz++) { - array(ispec, igllx, iglly, igllz) = dummy_T[counter]; + try { + // Read into dummy vector + for (int i = 0; i < n0; i++) { + specfem::IO::fortran_read_line(stream, &dummy_T); + // Assign to KokkosView + int counter = 0; + for (int j = 0; j < n1; j++) { + for (int k = 0; k < n2; k++) { + array(i, j, k) = dummy_T[counter]; counter++; } } } - if ((ispec == 0) || (ispec == nspec - 1)) { - std::cout << "ispec=" << ispec - << " array(ispec,0,0,0)=" << array(ispec, 0, 0, 0) - << std::endl; - } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 3D array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); } - - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading array from database file:\n" - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); } -} - -template -void specfem::IO::mesh::impl::fortran::dim3::read_array(std::ifstream &stream, - View5D &array) { - const int nspec = array.extent(0); - const int ngllx = array.extent(1); - const int nglly = array.extent(2); - const int ngllz = array.extent(3); - const int ncomp = array.extent(4); + // 4D array implementation + else if constexpr (rank == 4) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); + const int n2 = array.extent(2); + const int n3 = array.extent(3); - std::vector dummy_T(ngllx * nglly * ngllz * ncomp, -9999999); - - try { - for (int ispec = 0; ispec < nspec; ispec++) { - // Read into dummy vector for each ispec - specfem::IO::fortran_read_line(stream, &dummy_T); + std::vector dummy_T(n1 * n2 * n3, -99999); - // Assign to KokkosView - int counter = 0; - for (int igllx = 0; igllx < ngllx; igllx++) { - for (int iglly = 0; iglly < nglly; iglly++) { - for (int igllz = 0; igllz < ngllz; igllz++) { - for (int icomp = 0; icomp < ncomp; icomp++) { - array(ispec, igllx, iglly, igllz, icomp) = dummy_T[counter]; + try { + for (int i = 0; i < n0; i++) { + // Read into dummy vector + specfem::IO::fortran_read_line(stream, &dummy_T); + // Assign to KokkosView + int counter = 0; + for (int j = 0; j < n1; j++) { + for (int k = 0; k < n2; k++) { + for (int l = 0; l < n3; l++) { + array(i, j, k, l) = dummy_T[counter]; counter++; } } } } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 4D array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); } - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading array from database file:\n" - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); } -} + // 5D array implementation + else if constexpr (rank == 5) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); + const int n2 = array.extent(2); + const int n3 = array.extent(3); + const int n4 = array.extent(4); -template -void specfem::IO::mesh::impl::fortran::dim3::read_index_array( - std::ifstream &stream, View1D &array) { - const int n = array.extent(0); - std::vector dummy_T(n, -999999); + std::vector dummy_T(n1 * n2 * n3 * n4, -99999); - try { - // Read into dummy vector - specfem::IO::fortran_read_line(stream, &dummy_T); - - // Assign to KokkosView - for (int i = 0; i < n; i++) { - array(i) = dummy_T[i] - 1; + try { + for (int i = 0; i < n0; i++) { + // Read into dummy vector + specfem::IO::fortran_read_line(stream, &dummy_T); + // Assign to KokkosView + int counter = 0; + for (int j = 0; j < n1; j++) { + for (int k = 0; k < n2; k++) { + for (int l = 0; l < n3; l++) { + for (int m = 0; m < n4; m++) { + array(i, j, k, l, m) = dummy_T[counter]; + counter++; + } + } + } + } + } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 5D array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); } - - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading 1D index_array from database file:\n" - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); + } else { + throw std::runtime_error("Unsupported rank for read_array array"); } -} +}; -template +template void specfem::IO::mesh::impl::fortran::dim3::read_index_array( - std::ifstream &stream, View4D &array) { - const int nspec = array.extent(0); - const int ngllx = array.extent(1); - const int nglly = array.extent(2); - const int ngllz = array.extent(3); + std::ifstream &stream, ViewType &array) { - std::vector dummy_T(ngllx * nglly * ngllz, -9999.0); + // Get value_type and rank of the ViewType + using value_type = typename ViewType::value_type; + constexpr int rank = ViewType::rank; - try { - for (int ispec = 0; ispec < nspec; ispec++) { - // Read into dummy vector for each ispec + // 1D array implementation + if constexpr (rank == 1) { + + const int n = array.extent(0); + std::vector dummy_T(n, -999999); + + try { + // Read into dummy vector specfem::IO::fortran_read_line(stream, &dummy_T); // Assign to KokkosView - int counter = 0; - for (int igllx = 0; igllx < ngllx; igllx++) { - for (int iglly = 0; iglly < nglly; iglly++) { - for (int igllz = 0; igllz < ngllz; igllz++) { - array(ispec, igllx, iglly, igllz) = dummy_T[counter] - 1; + for (int i = 0; i < n; i++) { + array(i) = dummy_T[i] - 1; + } + + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 1D index_array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } + } + // 2D array implementation + else if constexpr (rank == 2) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); + + std::vector dummy_T(n1, -999999); + + try { + for (int i = 0; i < n0; i++) { + // Read into dummy vector + specfem::IO::fortran_read_line(stream, &dummy_T); + + // Assign to KokkosView + int counter = 0; + for (int j = 0; j < n1; j++) { + array(i, j) = dummy_T[counter] - 1; + counter++; + } + } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 2D index_array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } + } + // 3D array implementation + else if constexpr (rank == 3) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); + const int n2 = array.extent(2); + + std::vector dummy_T(n1 * n2, -999999); + + try { + for (int i = 0; i < n0; i++) { + // Read into dummy vector + specfem::IO::fortran_read_line(stream, &dummy_T); + + // Assign to KokkosView + int counter = 0; + for (int j = 0; j < n1; j++) { + for (int k = 0; k < n2; k++) { + array(i, j, k) = dummy_T[counter] - 1; counter++; } } } - if ((ispec == 0) || (ispec == nspec - 1)) { - std::cout << "ispec=" << ispec - << " array(ispec,0,0,0)=" << array(ispec, 0, 0, 0) - << std::endl; + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 3D index_array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } + } + // 4D array implementation + else if constexpr (rank == 4) { + const int n0 = array.extent(0); + const int n1 = array.extent(1); + const int n2 = array.extent(2); + const int n3 = array.extent(3); + + std::vector dummy_T(n1 * n2 * n3, -999999); + + try { + for (int i = 0; i < n0; i++) { + // Read into dummy vector + specfem::IO::fortran_read_line(stream, &dummy_T); + + // Assign to KokkosView + int counter = 0; + for (int j = 0; j < n1; j++) { + for (int k = 0; k < n2; k++) { + for (int l = 0; l < n3; l++) { + array(i, j, k, l) = dummy_T[counter] - 1; + counter++; + } + } + } } + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading 4D index_array from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); } + } else { + throw std::runtime_error("Unsupported rank for read_index_array"); + } +}; + +namespace specfem { +namespace IO { +namespace mesh { +namespace impl { +namespace fortran { +namespace dim3 { - } catch (std::runtime_error &e) { +// Primary template for try-catch wrapper +/** + * @brief Try-catch wrapper for fortran_read_line + * + * @param message Message to print if exception is caught + * @param args Arguments to pass to fortran_read_line + * @throws std::runtime_error if an error occurs while reading the line + * includes the input message, so the user know which value errors + * @see ::specfem::IO::fortran_read_line + * @note This function is a wrapper for fortran_read_line that catches + * exceptions and throws a runtime_error with a more informative message + * + * @code{.cpp} + * // Example of how to use this function + * int very_specific_value; + * try_read_line("very_specific_value", stream, &value); + * @endcode + */ +template +auto try_read_line(const std::string &message, Args &&...args) + -> decltype(specfem::IO::fortran_read_line(std::forward(args)...)) { + try { + return specfem::IO::fortran_read_line(std::forward(args)...); + } catch (const std::exception &e) { + std::ostringstream error_message; + error_message << "Exception caught in try_read_line(" << message + << ", ...): " << e.what(); + throw std::runtime_error(error_message.str()); + } catch (...) { std::ostringstream error_message; - error_message << "Error reading array from database file:\n" - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + error_message << "Unknown exception caught in " << message; throw std::runtime_error(error_message.str()); } } + +/** + * @brief Try-catch wrapper for read_array + * + * @param message Message to print if exception is caught + * @param args Arguments to pass to read_array + * @throws std::runtime_error if an error occurs while reading the array + * includes the input message, so the user know which array errors + * @see ::specfem::IO::mesh::impl::fortran::dim3::read_array + * @note This function is a wrapper for read_array that catches exceptions + * and throws a runtime_error with a more informative message + * + * @code{.cpp} + * // Example of how to use this function + * Kokkos::View specific_array("array", 10); + * try_read_array("specific_array", stream, array); + * @endcode + */ +template +auto try_read_array(const std::string &message, Args &&...args) + -> decltype(specfem::IO::mesh::impl::fortran::dim3::read_array( + std::forward(args)...)) { + try { + return specfem::IO::mesh::impl::fortran::dim3::read_array( + std::forward(args)...); + } catch (const std::exception &e) { + std::ostringstream error_message; + error_message << "Exception caught in try_read_array('" << message + << "', ...): " << e.what(); + throw std::runtime_error(error_message.str()); + } catch (...) { + std::ostringstream error_message; + error_message << "Unknown exception caught in try_read_array('" << message + << "', ...)."; + throw std::runtime_error(error_message.str()); + } +} + +/** + * @brief Try-catch wrapper for read_index_array + * + * @param message Message to print if exception is caught + * @param args Arguments to pass to read_index_array + * @throws std::runtime_error if an error occurs while reading the array that + * includes the input message, so the user know which array errors + * @see ::specfem::IO::mesh::impl::fortran::dim3::read_index_array + * @note This function is a wrapper for read_index_array that catches exceptions + * and throws a runtime_error with a more informative message + * + * @code{.cpp} + * // Example of how to use this function + * Kokkos::View specific_array("array", 10); + * try_read_index_array("specific_array", stream, array); + * @endcode + */ +template +auto try_read_index_array(const std::string &message, Args &&...args) + -> decltype(specfem::IO::mesh::impl::fortran::dim3::read_index_array( + std::forward(args)...)) { + try { + return specfem::IO::mesh::impl::fortran::dim3::read_index_array( + std::forward(args)...); + } catch (const std::exception &e) { + std::ostringstream error_message; + error_message << "Exception caught in try_read_index('" << message + << "', ...): " << e.what(); + throw std::runtime_error(error_message.str()); + } catch (...) { + std::ostringstream error_message; + error_message << "Unknown exception caught in try_read_index('" << message + << "', ...)."; + throw std::runtime_error(error_message.str()); + } +} + +} // namespace dim3 +} // namespace fortran +} // namespace impl +} // namespace mesh +} // namespace IO +} // namespace specfem diff --git a/include/compute/fields/data_access.tpp b/include/compute/fields/data_access.tpp index 1fc1cae91..94f27af49 100644 --- a/include/compute/fields/data_access.tpp +++ b/include/compute/fields/data_access.tpp @@ -20,17 +20,7 @@ KOKKOS_FORCEINLINE_FUNCTION void impl_load_on_device(const int iglob, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < ViewType::components; ++icomp) { @@ -79,17 +69,7 @@ impl_load_on_device(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { @@ -152,17 +132,7 @@ impl_load_on_device(const specfem::point::simd_assembly_index &index, mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < components; ++icomp) { @@ -208,17 +178,7 @@ inline void impl_load_on_host(const int iglob, const WavefieldType &field, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < ViewType::components; ++icomp) { @@ -266,17 +226,7 @@ inline void impl_load_on_host(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { @@ -338,17 +288,7 @@ inline void impl_load_on_host(const specfem::point::simd_assembly_index &index, mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < components; ++icomp) { @@ -501,17 +441,7 @@ impl_store_on_device(const int iglob, const ViewType &point_field, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); + const auto &curr_field = field.template get_field(); for (int icomp = 0; icomp < ViewType::components; ++icomp) { if constexpr (StoreDisplacement) { @@ -547,18 +477,7 @@ impl_store_on_device(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { continue; @@ -620,22 +539,11 @@ impl_store_on_device(const specfem::point::simd_assembly_index &index, mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < components; ++icomp) { Kokkos::Experimental::where(mask, point_field.displacement(icomp)) - .copy_to(&curr_field.h_field(iglob, icomp), tag_type()); + .copy_to(&curr_field.field(iglob, icomp), tag_type()); } } @@ -676,18 +584,7 @@ inline void impl_store_on_host(const int iglob, const ViewType &point_field, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int icomp = 0; icomp < ViewType::components; ++icomp) { if constexpr (StoreDisplacement) { curr_field.h_field(iglob, icomp) = point_field.displacement(icomp); @@ -722,18 +619,7 @@ inline void impl_store_on_host(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { continue; @@ -795,18 +681,7 @@ inline void impl_store_on_host(const specfem::point::simd_assembly_index &index, mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < components; ++icomp) { Kokkos::Experimental::where(mask, point_field.displacement(icomp)) @@ -968,18 +843,7 @@ impl_add_on_device(const int iglob, const ViewType &point_field, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int icomp = 0; icomp < ViewType::components; ++icomp) { if constexpr (StoreDisplacement) { curr_field.field(iglob, icomp) += point_field.displacement(icomp); @@ -1014,18 +878,7 @@ impl_add_on_device(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { continue; @@ -1087,18 +940,7 @@ impl_add_on_device(const specfem::point::simd_assembly_index &index, mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < components; ++icomp) { typename ViewType::simd::datatype lhs; @@ -1163,18 +1005,7 @@ inline void impl_add_on_host(const int iglob, const ViewType &point_field, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int icomp = 0; icomp < ViewType::components; ++icomp) { if constexpr (StoreDisplacement) { curr_field.h_field(iglob, icomp) += point_field.displacement(icomp); @@ -1209,18 +1040,7 @@ inline void impl_add_on_host(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { continue; @@ -1281,18 +1101,7 @@ inline void impl_add_on_host(const specfem::point::simd_assembly_index &index, mask_type mask([&](std::size_t lane) { return index.mask(lane); }); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); if constexpr (StoreDisplacement) { for (int icomp = 0; icomp < components; ++icomp) { typename ViewType::simd::datatype lhs; @@ -1472,18 +1281,7 @@ impl_atomic_add_on_device(const int iglob, const ViewType &point_field, constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int icomp = 0; icomp < ViewType::components; ++icomp) { if constexpr (StoreDisplacement) { Kokkos::atomic_add(&curr_field.field(iglob, icomp), @@ -1522,18 +1320,7 @@ impl_atomic_add_on_device(const typename ViewType::simd::mask_type &mask, constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { continue; @@ -1586,18 +1373,7 @@ inline void impl_atomic_add_on_host(const int iglob, const ViewType &point_field constexpr static bool StoreMassMatrix = ViewType::store_mass_matrix; constexpr static auto MediumType = ViewType::medium_tag; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int icomp = 0; icomp < ViewType::components; ++icomp) { if constexpr (StoreDisplacement) { Kokkos::atomic_add(&curr_field.h_field(iglob, icomp), @@ -1635,18 +1411,7 @@ inline void impl_atomic_add_on_host(const typename ViewType::simd::mask_type &ma constexpr static auto MediumType = ViewType::medium_tag; constexpr static int components = ViewType::components; - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); for (int lane = 0; lane < ViewType::simd::size(); ++lane) { if (!mask[lane]) { continue; @@ -1793,18 +1558,7 @@ impl_load_on_device(const MemberType &team, const int ispec, Kokkos::DefaultExecutionSpace>, "Calling team must have a device execution space"); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); Kokkos::parallel_for( Kokkos::TeamThreadRange(team, NGLL * NGLL), [&](const int &xz) { int iz, ix; @@ -1855,18 +1609,7 @@ inline void impl_load_on_host(const MemberType &team, const int ispec, std::is_same_v, "Calling team must have a host execution space"); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); Kokkos::parallel_for( Kokkos::TeamThreadRange(team, NGLL * NGLL), [&](const int &xz) { int iz, ix; @@ -1929,18 +1672,7 @@ impl_load_on_device(const MemberType &team, const IteratorType &iterator, typename ViewType::memory_space>::accessible, "Calling team must have access to the memory space of the view"); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); Kokkos::parallel_for( Kokkos::TeamThreadRange(team, iterator.chunk_size()), [&](const int &i) { const auto iterator_index = iterator(i); @@ -2006,18 +1738,7 @@ impl_load_on_device(const MemberType &team, const IteratorType &iterator, Kokkos::DefaultExecutionSpace>, "Calling team must have a device execution space"); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); Kokkos::parallel_for( Kokkos::TeamThreadRange(team, iterator.chunk_size()), [&](const int &i) { const auto iterator_index = iterator(i); @@ -2089,18 +1810,7 @@ inline void impl_load_on_host(const MemberType &team, const IteratorType &iterat typename ViewType::memory_space>::accessible, "Calling team must have access to the memory space of the view"); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); Kokkos::parallel_for( Kokkos::TeamThreadRange(team, iterator.chunk_size()), [&](const int &i) { const auto iterator_index = iterator(i); @@ -2164,18 +1874,7 @@ inline void impl_load_on_host(const MemberType &team, const IteratorType &iterat typename ViewType::memory_space>::accessible, "Calling team must have access to the memory space of the view"); - const auto &curr_field = - [&]() -> const specfem::compute::impl::field_impl< - specfem::dimension::type::dim2, MediumType> & { - if constexpr (MediumType == specfem::element::medium_tag::elastic) { - return field.elastic; - } else if constexpr (MediumType == specfem::element::medium_tag::acoustic) { - return field.acoustic; - } else { - static_assert("medium type not supported"); - } - }(); - + const auto &curr_field = field.template get_field(); Kokkos::parallel_for( Kokkos::TeamThreadRange(team, iterator.chunk_size()), [&](const int &i) { const auto iterator_index = iterator(i); diff --git a/include/compute/fields/simulation_field.hpp b/include/compute/fields/simulation_field.hpp index 70335c809..2fe7fb9b2 100644 --- a/include/compute/fields/simulation_field.hpp +++ b/include/compute/fields/simulation_field.hpp @@ -94,6 +94,23 @@ struct simulation_field { } } + /** + * @brief Returns the field for a given medium + * + */ + template + KOKKOS_INLINE_FUNCTION constexpr specfem::compute::impl::field_impl< + specfem::dimension::type::dim2, MediumTag> const & + get_field() const { + if constexpr (MediumTag == specfem::element::medium_tag::elastic) { + return elastic; + } else if constexpr (MediumTag == specfem::element::medium_tag::acoustic) { + return acoustic; + } else { + static_assert("medium type not supported"); + } + } + int nglob = 0; ///< Number of global degrees of freedom int nspec; ///< Number of spectral elements int ngllz; ///< Number of quadrature points in z direction diff --git a/include/compute/properties/properties.hpp b/include/compute/properties/properties.hpp index 2addf56aa..d0ba31bf0 100644 --- a/include/compute/properties/properties.hpp +++ b/include/compute/properties/properties.hpp @@ -2,6 +2,7 @@ #include "compute/element_types/element_types.hpp" #include "compute/impl/value_containers.hpp" +#include "enumerations/dimension.hpp" #include "enumerations/specfem_enums.hpp" #include "kokkos_abstractions.h" #include "macros.hpp" @@ -47,9 +48,11 @@ struct properties * @param has_gll_model Whether a GLL model is present (skip material property * assignment if true) */ - properties(const int nspec, const int ngllz, const int ngllx, - const specfem::compute::element_types &element_types, - const specfem::mesh::materials &materials, bool has_gll_model); + properties( + const int nspec, const int ngllz, const int ngllx, + const specfem::compute::element_types &element_types, + const specfem::mesh::materials &materials, + bool has_gll_model); ///@} diff --git a/include/medium/material_properties.hpp b/include/medium/material_properties.hpp index cc8d7df1a..46e9f267b 100644 --- a/include/medium/material_properties.hpp +++ b/include/medium/material_properties.hpp @@ -21,7 +21,8 @@ struct material_properties material_properties( const Kokkos::View elements, const int ngllz, const int ngllx, - const specfem::mesh::materials &materials, const bool has_gll_model, + const specfem::mesh::materials &materials, + const bool has_gll_model, const specfem::kokkos::HostView1d property_index_mapping) : specfem::medium::properties_container( elements.extent(0), ngllz, ngllx) { diff --git a/include/mesh/dim2/coupled_interfaces/coupled_interfaces.hpp b/include/mesh/dim2/coupled_interfaces/coupled_interfaces.hpp index b3508ece5..4d825fd8b 100644 --- a/include/mesh/dim2/coupled_interfaces/coupled_interfaces.hpp +++ b/include/mesh/dim2/coupled_interfaces/coupled_interfaces.hpp @@ -1,21 +1,24 @@ #pragma once +#include "enumerations/dimension.hpp" #include "enumerations/medium.hpp" #include "interface_container.hpp" +#include "mesh/mesh_base.hpp" #include "specfem_mpi/specfem_mpi.hpp" #include #include namespace specfem { namespace mesh { + /** - * @brief Information about coupled interfaces + * @brief Struct to store coupled interfaces for the 2D mesh * */ -template -struct coupled_interfaces { +template <> struct coupled_interfaces { public: - constexpr static auto dimension = DimensionType; ///< Dimension + constexpr static auto dimension = + specfem::dimension::type::dim2; ///< Dimension /** * @name Constructors * @@ -29,15 +32,15 @@ struct coupled_interfaces { : elastic_acoustic(), acoustic_poroelastic(), elastic_poroelastic(){}; coupled_interfaces(specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::elastic, + dimension, specfem::element::medium_tag::elastic, specfem::element::medium_tag::acoustic> elastic_acoustic, specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::acoustic, + dimension, specfem::element::medium_tag::acoustic, specfem::element::medium_tag::poroelastic> acoustic_poroelastic, specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::elastic, + dimension, specfem::element::medium_tag::elastic, specfem::element::medium_tag::poroelastic> elastic_poroelastic) : elastic_acoustic(elastic_acoustic), @@ -55,27 +58,27 @@ struct coupled_interfaces { template std::variant, specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::acoustic, + dimension, specfem::element::medium_tag::acoustic, specfem::element::medium_tag::poroelastic>, specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::elastic, + dimension, specfem::element::medium_tag::elastic, specfem::element::medium_tag::poroelastic> > get() const; - specfem::mesh::interface_container elastic_acoustic; ///< Elastic-acoustic interfaces - specfem::mesh::interface_container acoustic_poroelastic; ///< Acoustic-poroelastic interfaces - specfem::mesh::interface_container elastic_poroelastic; ///< Elastic-poroelastic interfaces diff --git a/include/mesh/dim2/materials/materials.hpp b/include/mesh/dim2/materials/materials.hpp index fec3b49a5..40c5adc5c 100644 --- a/include/mesh/dim2/materials/materials.hpp +++ b/include/mesh/dim2/materials/materials.hpp @@ -1,7 +1,9 @@ #pragma once +#include "enumerations/dimension.hpp" #include "kokkos_abstractions.h" #include "medium/material.hpp" +#include "mesh/mesh_base.hpp" #include "specfem_mpi/interface.hpp" #include "specfem_setup.hpp" #include @@ -12,7 +14,10 @@ namespace mesh { * @brief Material properties information * */ -struct materials { + +template <> struct materials { + constexpr static auto dimension = + specfem::dimension::type::dim2; ///< Dimension type struct material_specification { specfem::element::medium_tag type; ///< Type of element @@ -59,17 +64,19 @@ struct materials { material_index_mapping; ///< Mapping of spectral element to material ///< properties - specfem::mesh::materials::material + specfem::mesh::materials::material< + specfem::element::medium_tag::elastic, + specfem::element::property_tag::isotropic> elastic_isotropic; ///< Elastic isotropic material properties - specfem::mesh::materials::material< + specfem::mesh::materials::material< specfem::element::medium_tag::elastic, specfem::element::property_tag::anisotropic> elastic_anisotropic; ///< Elastic anisotropic material properties - specfem::mesh::materials::material + specfem::mesh::materials::material< + specfem::element::medium_tag::acoustic, + specfem::element::property_tag::isotropic> acoustic_isotropic; ///< Acoustic isotropic material properties /** diff --git a/include/mesh/dim2/materials/materials.tpp b/include/mesh/dim2/materials/materials.tpp index 867ed0389..42653d5e7 100644 --- a/include/mesh/dim2/materials/materials.tpp +++ b/include/mesh/dim2/materials/materials.tpp @@ -9,7 +9,7 @@ template -specfem::mesh::materials::material::material( + specfem::mesh::materials::material::material( const int n_materials, const std::vector > &l_materials) diff --git a/include/mesh/dim2/mesh.hpp b/include/mesh/dim2/mesh.hpp index 5ec58959b..d4c099c9b 100644 --- a/include/mesh/dim2/mesh.hpp +++ b/include/mesh/dim2/mesh.hpp @@ -60,7 +60,8 @@ template <> struct mesh { ///< nodes ///< (never ///< used) - specfem::mesh::materials materials; ///< Defines material properties + specfem::mesh::materials materials; ///< Defines material + ///< properties /** * @name Constructors @@ -73,27 +74,69 @@ template <> struct mesh { */ mesh(){}; - mesh(const int npgeo, const int nspec, const int nproc, - const specfem::mesh::control_nodes - &control_nodes, - const specfem::mesh::parameters - ¶meters, - const specfem::mesh::coupled_interfaces - &coupled_interfaces, - const specfem::mesh::boundaries - &boundaries, - const specfem::mesh::tags &tags, - const specfem::mesh::elements::tangential_elements< - specfem::dimension::type::dim2> &tangential_nodes, - const specfem::mesh::elements::axial_elements< - specfem::dimension::type::dim2> &axial_nodes, - const specfem::mesh::materials &materials) + /** + * @brief Mesh constructor + * + * This constructor initializes the mesh struct with the given parameters + * + * @param npgeo Total number of spectral element control nodes + * @param nspec Total number of spectral elements + * @param nproc Total number of processors + * @param control_nodes Struct to store control nodes + * @param parameters Struct to store simulation launch parameters + * @param coupled_interfaces Struct to store coupled interfaces + * @param boundaries Struct to store information at the boundaries + * @param tags Struct to store tags for every spectral element + * @param tangential_nodes Struct to store tangential nodes + * @param axial_nodes Struct to store axial nodes + * @param materials Struct to store material properties + * + * @see ::specfem::mesh::control_nodes, ::specfem::mesh::parameters, + * ::specfem::mesh::coupled_interfaces, ::specfem::mesh::boundaries, + * ::specfem::mesh::tags, ::specfem::mesh::elements::tangential_elements, + * ::specfem::mesh::elements::axial_elements, ::specfem::mesh::materials + * + * @code{.cpp} + * // Example of how to use this constructor + * specfem::mesh::mesh mesh( + * npgeo, nspec, nproc, control_nodes, parameters, coupled_interfaces, + * boundaries, tags, tangential_nodes, axial_nodes, materials); + * @endcode + */ + mesh( + const int npgeo, const int nspec, const int nproc, + const specfem::mesh::control_nodes + &control_nodes, + const specfem::mesh::parameters + ¶meters, + const specfem::mesh::coupled_interfaces + &coupled_interfaces, + const specfem::mesh::boundaries + &boundaries, + const specfem::mesh::tags &tags, + const specfem::mesh::elements::tangential_elements< + specfem::dimension::type::dim2> &tangential_nodes, + const specfem::mesh::elements::axial_elements< + specfem::dimension::type::dim2> &axial_nodes, + const specfem::mesh::materials &materials) : npgeo(npgeo), nspec(nspec), nproc(nproc), control_nodes(control_nodes), parameters(parameters), coupled_interfaces(coupled_interfaces), boundaries(boundaries), tags(tags), tangential_nodes(tangential_nodes), axial_nodes(axial_nodes), materials(materials){}; - ///@} + ///@} // Constructors + /** + * @brief Print mesh information + * + * This function prints the mesh information + * + * @return std::string String containing the mesh information + * + * @code{.cpp} + * // Example of how to use this function + * std::string mesh_info = mesh.print(); + * @endcode + */ std::string print() const; }; } // namespace mesh diff --git a/include/mesh/dim2/tags/tags.hpp b/include/mesh/dim2/tags/tags.hpp index 78df66968..07844bfe5 100644 --- a/include/mesh/dim2/tags/tags.hpp +++ b/include/mesh/dim2/tags/tags.hpp @@ -41,9 +41,10 @@ template <> struct tags { * @param materials Material properties * @param boundaries Boundary information */ - tags(const specfem::mesh::materials &materials, - const specfem::mesh::boundaries - &boundaries); + tags( + const specfem::mesh::materials &materials, + const specfem::mesh::boundaries + &boundaries); ///@} }; } // namespace mesh diff --git a/include/mesh/dim3/boundaries/absorbing_boundary.hpp b/include/mesh/dim3/boundaries/absorbing_boundary.hpp new file mode 100644 index 000000000..62799b4c2 --- /dev/null +++ b/include/mesh/dim3/boundaries/absorbing_boundary.hpp @@ -0,0 +1,277 @@ +#pragma once + +#include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" +#include "mesh/mesh_base.hpp" +#include "specfem_setup.hpp" +#include + +namespace specfem { +namespace mesh { + +/** + * @brief Struct to store absorbing boundaries + * + */ +template struct stacey_mass {}; + +/** + * @brief Struct to store absorbing boundaries for an elastic medium + * + */ +template <> struct stacey_mass { + int nglob; ///< Number of GLL points + bool acoustic = false; ///< Flag for acoustic simulation + bool elastic = true; ///< Flag for elastic simulation + + Kokkos::View x; + Kokkos::View y; + Kokkos::View z; + + /** @name Stacey Elastic Mass struct constructors + * @{ + */ + /** + * @brief Default constructor + * + */ + stacey_mass(){}; + + /** + * @brief Constructor + * + * @param nglob Number of GLL points + * + * @note The x, y, and z views are created + * + * @code{.cpp} + * // Example of how to use this constructor + * const int nglob = 10; + * stacey_mass mass_matrix(nglob); + * + * // Populate the mass matrix from the binary file + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, mass_matrix.x); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, mass_matrix.y); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, mass_matrix.z); + * @endcode + */ + stacey_mass(const int nglob) : nglob(nglob) { + x = Kokkos::View("rmass_x", nglob); + y = Kokkos::View("rmass_y", nglob); + z = Kokkos::View("rmass_z", nglob); + } + /** @} */ // Stacey Elastic Mass struct constructors +}; + +/** + * @brief Struct to store absorbing boundaries for an acoustic medium + * + */ +template <> struct stacey_mass { + int nglob; ///< Number of GLL points + bool acoustic = true; ///< Flag for acoustic simulation + bool elastic = false; ///< Flag for elastic simulation + + Kokkos::View mass; + + /** @name Constructors + * @{ + */ + + /** + * @brief Default constructor + * + */ + stacey_mass(){}; + + /** + * @brief Constructor + * + * @param nglob Number of GLL points + * + * @note The mass view is created with the name "rmass_x" + * + * @code{.cpp} + * // Example of how to use this constructor + * const int nglob = 10; + * stacey_mass mass_matrix(nglob); + * + * // Populate the mass matrix from the binary file + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, + * mass_matrix.mass); + * @endcode + */ + stacey_mass(const int nglob) : nglob(nglob) { + mass = Kokkos::View("rmass_x", nglob); + } + /** @} */ // Constructors +}; + +/** + * @brief Struct to store absorbing boundaries for 3D meshes + */ +template <> struct absorbing_boundary { + + constexpr static auto dimension = + specfem::dimension::type::dim3; ///< Dimension type + // const int ndim = 3; ///< Dimension + + int num_abs_boundary_faces; ///< Number of boundaries faces + int nelements; ///< Number of elements on the boundary + int ngllsquare; ///< Number of GLL points squared + bool acoustic = false; ///< Flag for acoustic simulation + bool elastic = false; ///< Flag for elastic simulation + + int nspec2D_xmin, nspec2D_xmax; /// Number of elements on the boundaries + int nspec2D_ymin, nspec2D_ymax; + int NSPEC2D_BOTTOM, NSPEC2D_TOP; + + Kokkos::View ibelm_xmin, + ibelm_xmax; ///< Spectral + ///< element + ///< index for + ///< elements on + ///< the boundary + Kokkos::View ibelm_ymin, ibelm_ymax; + Kokkos::View ibelm_bottom, ibelm_top; + + Kokkos::View ispec; ///< Spectral element index for + ///< elements on the boundary + Kokkos::View ijk; ///< Which edge of the element + ///< is on the boundary + Kokkos::View jacobian2Dw; ///< Jacobian of + ///< the 2D + Kokkos::View normal; ///< Jacobian of the 2D + + stacey_mass mass_elastic; + stacey_mass mass_acoustic; + + /** @name Constructors + * @{ + */ + + /** + * @brief Default constructor + * + */ + absorbing_boundary(){}; + + /** + * @brief Constructor for absorbing boundaries + * + * This struct holds views with the names "ispec", "ijk", "jacobian2Dw", + * and "normal" and the mass matrices are created if the medium is elastic or + * acoustic under the names "mass_elastic" and "mass_acoustic" respectively. + * + * @param nglob Number of GLL points + * @param num_abs_boundary_faces Number of boundary faces + * @param ngllsquare Number of GLL points squared + * @param acoustic Flag for acoustic simulation + * @param elastic Flag for elastic simulation + * @param nspec2D_xmin Number of elements on the x-min boundary + * @param nspec2D_xmax Number of elements on the x-max boundary + * @param nspec2D_ymin Number of elements on the y-min boundary + * @param nspec2D_ymax Number of elements on the y-max boundary + * @param NSPEC2D_BOTTOM Number of elements on the bottom boundary + * @param NSPEC2D_TOP Number of elements on the top boundary + * + * + * @code{.cpp} + * // Example of how to use this constructor + * const int nglob = 10; + * const int num_abs_boundary_faces = 5; + * const int ngllsquare = 100; + * const bool acoustic = true; + * const bool elastic = false; + * const int nspec2D_xmin = 10; + * const int nspec2D_xmax = 10; + * const int nspec2D_ymin = 10; + * const int nspec2D_ymax = 10; + * const int NSPEC2D_BOTTOM = 10; + * const int NSPEC2D_TOP = 10; + * + * absorbing_boundary abs_boundary(nglob, + * num_abs_boundary_faces, ngllsquare, acoustic, elastic, nspec2D_xmin, + * nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, + * NSPEC2D_BOTTOM, NSPEC2D_TOP); + * + * // Populate the views from the binary file + * specfem::IO::mesh::impl::fortran::dim3::read_index_array(stream, + * abs_boundary.ispec); + * specfem::IO::mesh::impl::fortran::dim3::read_index_array(stream, + * abs_boundary.ijk); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, + * abs_boundary.jacobian2Dw); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, + * abs_boundary.normal); + * @endcode + */ + absorbing_boundary(const int nglob, const int num_abs_boundary_faces, + const int ngllsquare, const bool acoustic, + const bool elastic, const int nspec2D_xmin, + const int nspec2D_xmax, const int nspec2D_ymin, + const int nspec2D_ymax, const int NSPEC2D_BOTTOM, + const int NSPEC2D_TOP) + : nelements(num_abs_boundary_faces), ngllsquare(ngllsquare), + num_abs_boundary_faces(num_abs_boundary_faces), acoustic(acoustic), + elastic(elastic), nspec2D_xmin(nspec2D_xmin), + nspec2D_xmax(nspec2D_xmax), nspec2D_ymin(nspec2D_ymin), + nspec2D_ymax(nspec2D_ymax), NSPEC2D_BOTTOM(NSPEC2D_BOTTOM), + NSPEC2D_TOP(NSPEC2D_TOP) { + + ispec = Kokkos::View("ispec", nelements); + ijk = Kokkos::View("ijk", nelements, 3, + ngllsquare); + jacobian2Dw = Kokkos::View( + "jacobian2Dw", nelements, ngllsquare); + // ndim=3 + + normal = Kokkos::View("normal", nelements, + 3, ngllsquare); + + if (elastic) { + mass_elastic = stacey_mass(nglob); + } + if (acoustic) { + mass_acoustic = + stacey_mass(nglob); + } + + // Create the boundary view elements only if the number of elements is + // greater than 0 on said boundary + if (nspec2D_xmin > 0) { + ibelm_xmin = + Kokkos::View("ibelm_xmin", nspec2D_xmin); + } + if (nspec2D_xmax > 0) { + ibelm_xmax = + Kokkos::View("ibelm_xmax", nspec2D_xmax); + } + if (nspec2D_ymin > 0) { + ibelm_ymin = + Kokkos::View("ibelm_ymin", nspec2D_ymin); + } + if (nspec2D_ymax > 0) { + ibelm_ymax = + Kokkos::View("ibelm_ymax", nspec2D_ymax); + } + if (NSPEC2D_BOTTOM > 0) { + ibelm_bottom = Kokkos::View("ibelm_bottom", + NSPEC2D_BOTTOM); + } + if (NSPEC2D_TOP > 0) { + ibelm_top = + Kokkos::View("ibelm_top", NSPEC2D_TOP); + } + } + /** @} */ // Constructors + + /** + * @brief Print basic information on the absorbing boundary struct + * + */ + std::string print() const; +}; + +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/boundaries/boundaries.hpp b/include/mesh/dim3/boundaries/boundaries.hpp new file mode 100644 index 000000000..26217551f --- /dev/null +++ b/include/mesh/dim3/boundaries/boundaries.hpp @@ -0,0 +1,14 @@ +#pragma once + +#include "enumerations/dimension.hpp" +#include "mesh/mesh_base.hpp" +#include "specfem_setup.hpp" +#include + +namespace specfem { +namespace mesh { + +template <> struct boundaries {}; + +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/boundaries/free_surface.hpp b/include/mesh/dim3/boundaries/free_surface.hpp new file mode 100644 index 000000000..51632fdee --- /dev/null +++ b/include/mesh/dim3/boundaries/free_surface.hpp @@ -0,0 +1,88 @@ +#pragma once + +#include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" +#include "mesh/mesh_base.hpp" +#include "specfem_setup.hpp" +#include + +namespace specfem { +namespace mesh { + +template <> struct free_surface { + + constexpr static auto dimension = + specfem::dimension::type::dim3; ///< Dimension type + // const int ndim = 3; ///< Dimension + + int num_free_surface_faces; ///< Number of boundaries faces + int nelements; ///< Number of elements on the boundary + int ngllsquare; ///< Number of GLL points squared + + Kokkos::View ispec; ///< Spectral element index for + ///< elements on the boundary + Kokkos::View ijk; ///< Which edge of the element + ///< is on the boundary + Kokkos::View jacobian2Dw; ///< Jacobian of + ///< the 2D + Kokkos::View normal; ///< Jacobian of the 2D + + /** + * @name Constructors + * + */ + ///@{ + /** + * @brief Default constructor + * + */ + free_surface(){}; + + /** + * @brief Constructor + * + * Constructor for the free_surface struct that initializes the number of + * the underlying arrays with the given parameters + * + * @param num_free_surface_faces Number of free surface faces + * @param ngllsquare Number of GLL points squared + * + * @code{.cpp} + * // Example of how to use this constructor + * specfem::mesh::free_surface free_surface( + * num_free_surface_faces, ngllsquare); + * @endcode + * + */ + free_surface(const int num_free_surface_faces, const int ngllsquare) + : nelements(num_free_surface_faces), ngllsquare(ngllsquare), + num_free_surface_faces(num_free_surface_faces) { + + ispec = Kokkos::View("ispec", nelements); + ijk = Kokkos::View("ijk", nelements, 3, + ngllsquare); + jacobian2Dw = Kokkos::View( + "jacobian2Dw", nelements, ngllsquare); + // ndim=3 + + normal = Kokkos::View("normal", nelements, + 3, ngllsquare); + } + ///@} + + /** + * @brief Print the free_surface struct + * + * Print the free_surface struct to the console + * + * @code{.cpp} + * // Example of how to use this function + * free_surface.print(); + * @endcode + * + */ + std::string print() const; +}; + +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/coordinates/coordinates.hpp b/include/mesh/dim3/coordinates/coordinates.hpp index ca72bc397..d78f9a476 100644 --- a/include/mesh/dim3/coordinates/coordinates.hpp +++ b/include/mesh/dim3/coordinates/coordinates.hpp @@ -1,13 +1,16 @@ #pragma once #include "enumerations/dimension.hpp" +#include "mesh/mesh_base.hpp" #include "specfem_setup.hpp" #include namespace specfem { namespace mesh { -template struct coordinates; - +/** + * @brief Struct to store coordinates for a 3D mesh + * + */ template <> struct coordinates { constexpr static auto dimension = specfem::dimension::type::dim3; @@ -25,19 +28,68 @@ template <> struct coordinates { UniqueView y; UniqueView z; - // Constructors + /** + * @name Constructors + * + * Constructors for the coordinates struct + * + */ + ///@{ + /** + * @brief Default constructor + * + */ coordinates(){}; // Default constructor - // Constructor to initialize the coordinates + /** + * @brief Constructor + * + * Constructor for the coordinates struct that initializes the number of + * the underlying arrays with the given parameters + * + * @param nspec Total number of spectral elements + * @param nglob Total number of global nodes + * @param ngllx Number of GLL points in the x direction + * @param nglly Number of GLL points in the y direction + * @param ngllz Number of GLL points in the z direction + * + * @code{.cpp} + * // Example of how to use this constructor + * specfem::mesh::coordinates coordinates( + * nspec, nglob, ngllx, nglly, ngllz); + * + * // Read coordinates from file + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, coordinates.x); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, coordinates.y); + * specfem::IO::mesh::impl::fortran::dim3::read_array(stream, coordinates.z); + * @endcode + * + */ coordinates(int nspec, int nglob, int ngllx, int nglly, int ngllz) : nspec(nspec), nglob(nglob), ngllx(ngllx), nglly(nglly), ngllz(ngllz), x("x", nglob), y("y", nglob), z("z", nglob){}; - void print() const; - - void print(int iglob) const; + ///@} // Constructors + /** + * @brief Print the coordinate specs + * + * @code{.cpp} + * // Example of how to use this function + * coordinates.print(); + * @endcode + */ + std::string print() const; - void print(int ispec, int igllx, int iglly, int igllz) const; + /** + * @brief Print the coordinates at a specific global node + * + * @code + * // Example of how to use this function + * int iglob = 10; + * coordinates.print(iglob); + * @endcode + */ + std::string print(int iglob) const; }; } // namespace mesh diff --git a/include/mesh/dim3/coupled_interfaces/coupled_interfaces.hpp b/include/mesh/dim3/coupled_interfaces/coupled_interfaces.hpp new file mode 100644 index 000000000..be02e6d9f --- /dev/null +++ b/include/mesh/dim3/coupled_interfaces/coupled_interfaces.hpp @@ -0,0 +1,191 @@ +#pragma once +#include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" +#include "mesh/mesh_base.hpp" + +namespace specfem { +namespace mesh { + +/** + * @brief Coupled interface struct that contains the information for the + * coupling between two different media + * + * @tparam MediumTag1 Medium tag for the first medium, e.g., elastic + * @tparam MediumTag2 Medium tag for the second medium, e.g., acoustic + * + */ +template +struct coupled_interface { + + constexpr static auto dimension = + specfem::dimension::type::dim3; ///< Dimension type + + int nelements; ///< Number of elements on the boundary + int ngllsquare; ///< Number of GLL points squared + int num_coupling_interface_faces; ///< Number of coupling interface faces + int nspec2D_xmin, nspec2D_xmax; ///< Number of elements on the boundaries + int nspec2D_ymin, nspec2D_ymax; ///< Number of elements on the boundaries + int NSPEC2D_BOTTOM, NSPEC2D_TOP; ///< Number of elements on the boundaries + + Kokkos::View ispec; ///< Spectral element index for + ///< elements on the boundary + Kokkos::View ijk; ///< Which edge of the element + ///< is on the boundary + Kokkos::View jacobian2Dw; ///< Jacobian of + ///< the 2D face + Kokkos::View normal; ///< Normal of the 2D + ///< face + + /** + * @name Constructors + * + */ + ///@{ + /** + * @brief Default constructor initializing an empty struct + * + */ + coupled_interface(){}; + + /** + * @brief Construct a new coupled interface object + * + * @param num_coupling_interface_faces + * @param ngllsquare + * + * @code{.cpp} + * // Example of how to use this constructor + * specfem::mesh::coupled_interface + * coupled_interface(num_coupling_interface_faces, ngllsquare); + * @endcode + */ + coupled_interface(const int num_coupling_interface_faces, + const int ngllsquare) + : nelements(num_coupling_interface_faces), ngllsquare(ngllsquare), + num_coupling_interface_faces(num_coupling_interface_faces) { + + ispec = Kokkos::View("ispec", nelements); + ijk = Kokkos::View("ijk", nelements, 3, + ngllsquare); + jacobian2Dw = Kokkos::View( + "jacobian2Dw", nelements, ngllsquare); + normal = Kokkos::View("normal", nelements, + 3, ngllsquare); + }; + + ///@} + + /** + * @brief Print basic information about the coupled interface + * + */ + void print() const; +}; + +/** + * @brief Coupled interfaces struct that contains the information for the + * coupling between multiple different media. + */ +template <> struct coupled_interfaces { + + constexpr static auto dimension = + specfem::dimension::type::dim3; ///< Dimension type + + bool acoustic_elastic = false; ///< Flag for acoustic elastic interfaces + bool acoustic_poroelastic = + false; ///< Flag for acoustic poroelastic interfaces + bool elastic_poroelastic = false; ///< Flag for elastic poroelastic interfaces + + int num_coupling_ac_el_faces; ///< Number of coupling acoustic elastic faces + int num_coupling_ac_po_faces; ///< Number of coupling acoustic poroelastic + ///< faces + int num_coupling_el_po_faces; ///< Number of coupling elastic poroelastic + ///< faces + int ngllsquare; ///< Number of GLL points squared + + // Create coupled placeholders. + coupled_interface + acoustic_elastic_interface; + coupled_interface + acoustic_poroelastic_interface; + coupled_interface + elastic_poroelastic_interface; + coupled_interface + poroelastic_elastic_interface; + + /** + * @name Constructors + * + */ + ///@{ + /** + * @brief Default constructor initializing an empty struct + * + */ + coupled_interfaces(){}; + + /** + * @brief Construct a new coupled interfaces object + * + + * + * @code{.cpp} + * // Example of how to use this constructor + * specfem::mesh::coupled_interfaces + * coupled_interfaces(num_coupling_ac_el_faces, num_coupling_ac_po_faces, + * num_coupling_el_po_faces, ngllsquare); + * @endcode + */ + coupled_interfaces(const int num_coupling_ac_el_faces, + const int num_coupling_ac_po_faces, + const int num_coupling_el_po_faces, const int ngllsquare) + : num_coupling_ac_el_faces(num_coupling_ac_el_faces), + num_coupling_ac_po_faces(num_coupling_ac_po_faces), + num_coupling_el_po_faces(num_coupling_el_po_faces), + ngllsquare(ngllsquare) { + + if (num_coupling_ac_el_faces > 0) { + acoustic_elastic_interface = + coupled_interface( + num_coupling_ac_el_faces, ngllsquare); + acoustic_elastic = true; + } + if (num_coupling_ac_po_faces > 0) { + acoustic_poroelastic_interface = + coupled_interface( + num_coupling_ac_po_faces, ngllsquare); + acoustic_poroelastic = true; + } + if (num_coupling_el_po_faces > 0) { + elastic_poroelastic_interface = + coupled_interface( + num_coupling_el_po_faces, ngllsquare); + + poroelastic_elastic_interface = + coupled_interface( + num_coupling_el_po_faces, ngllsquare); + + elastic_poroelastic = true; + } + }; + ///@} + + /** + * @brief Print basic information about the coupled interfaces to std::cout + * + */ + std::string print() const; +}; + +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/element_types/element_types.hpp b/include/mesh/dim3/element_types/element_types.hpp new file mode 100644 index 000000000..3594eddd9 --- /dev/null +++ b/include/mesh/dim3/element_types/element_types.hpp @@ -0,0 +1,118 @@ +#pragma once + +#include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" +#include "mesh/mesh_base.hpp" +#include +#include + +namespace specfem { +namespace mesh { + +/** + * @brief Struct to store element types for a 3D mesh + * + */ +template <> struct element_types { + + constexpr static auto dimension = specfem::dimension::type::dim3; + + int nspec; ///< Number of 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; + + /** + * @name Constructors + * + */ + ///@{ + /** + * @brief Default constructor initializing an empty struct + * + */ + element_types() = default; + + /** + * @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 + * int nspec = 10; + * specfem::mesh::element_types + * element_types(nspec); + * @endcode + */ + 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){}; + + ///@} + + /** + * @brief Set the elements for the given medium type + * + * This function loops over the boolean arrays ispec_is_elastic, + * ispec_is_acoustic, and ispec_is_poroelastic and creates arrays that contain + * the indices of the spectral elements that are elastic, acoustic, and + * poroelastic. + * + * @see get_elements + */ + void set_elements(View1D &ispec_type_in); + + /** + * @brief Get the elements for the given medium type + * + * @tparam MediumTag Medium type + * @return View1D Element indices for the given medium type + */ + template View1D get_elements(); + + /** + * @brief Print basic information about the element types + * + */ + std::string print() const; + + /** + * @brief Print the element type of a specific spectral element + * + * @param ispec Index of the spectral element + */ + std::string print(const int ispec) const; + + /** + * @brief Print the element type of element i/n + * + * @tparam MediumTag Medium type + */ + template + std::string print(const int i) const; + + 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 to global mapping + View1D ispec_acoustic; ///< Acoustic to global mapping + View1D ispec_poroelastic; ///< Poroelastic to global mapping +}; +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/mapping/mapping.hpp b/include/mesh/dim3/mapping/mapping.hpp index f73ab5d3f..335334e99 100644 --- a/include/mesh/dim3/mapping/mapping.hpp +++ b/include/mesh/dim3/mapping/mapping.hpp @@ -1,63 +1,87 @@ #pragma once #include "enumerations/dimension.hpp" +#include "mesh/mesh_base.hpp" #include "specfem_setup.hpp" #include namespace specfem { namespace mesh { -template struct mapping; - +/** + * @brief Struct to store mapping for a 3D mesh + * + */ template <> struct mapping { - constexpr static auto dimension = specfem::dimension::type::dim3; + constexpr static auto dimension = + specfem::dimension::type::dim3; ///< Dimension using UniqueViewInt = Kokkos::View; using UniqueViewBool = Kokkos::View; using LocalViewInt = Kokkos::View; // Parameters needed for ibool mapping - int nspec; - int nglob; - int nspec_irregular; + int nspec; ///< Number of spectral elements + int nglob; ///< Number of global elements + int nspec_irregular; ///< Number of irregular elements - int ngllx; - int nglly; - int ngllz; + int ngllx; ///< Number of GLL points in x direction + int nglly; ///< Number of GLL points in y direction + int ngllz; ///< Number of GLL points in z direction // I do not know currently what these are used for - type_real xix_regular; - type_real jacobian_regular; + type_real xix_regular; ///< Regular xi value + type_real jacobian_regular; ///< Regular Jacobian value // Indices of irregular elements size nspec_irregular - UniqueViewInt irregular_elements; + UniqueViewInt irregular_elements; ///< Irregular elements // ibool size nspec, ngllx, nglly, ngllz - LocalViewInt ibool; - - // Boolean array size nspec - UniqueViewBool is_acoustic; - UniqueViewBool is_elastic; - UniqueViewBool is_poroelastic; - - // Constructors + LocalViewInt ibool; ///< The local to global mapping + + /** + * @name Constructors + * + */ + ///@{ + /** + * @brief Construct a new mapping object + * + */ mapping(){}; // Default constructor - // Constructor to initialize the mapping + /** + * @brief Construct a new mapping object + * + * @param nspec Number of spectral elements + * @param nglob Number of global nodes + * @param nspec_irregular Number of irregular elements + * @param ngllx Number of GLL points in x direction + * @param nglly Number of GLL points in y direction + * @param ngllz Number of GLL points in z direction + */ mapping(int nspec, int nglob, int nspec_irregular, int ngllx, int nglly, int ngllz) : nspec(nspec), nglob(nglob), nspec_irregular(nspec_irregular), ngllx(ngllx), nglly(nglly), ngllz(ngllz), xix_regular(0.0), jacobian_regular(0.0), irregular_elements("irregular_elements", nspec_irregular), - ibool("ibool", nspec, ngllx, nglly, ngllz), - is_acoustic("is_acoustic", nspec), is_elastic("is_elastic", nspec), - is_poroelastic("is_poroelastic", nspec){}; - - void print() const; - - void print(int ispec) const; + ibool("ibool", nspec, ngllx, nglly, ngllz){}; + ///@} + + /** + * @brief Print basic information about the mapping + * + */ + std::string print() const; + + /** + * @brief Print the mapping for the given spectral element + * + * @param ispec Spectral element index + */ + std::string print(const int ispec) const; }; } // namespace mesh diff --git a/include/mesh/dim3/mass_matrix/mass_matrix.hpp b/include/mesh/dim3/mass_matrix/mass_matrix.hpp new file mode 100644 index 000000000..8ab548289 --- /dev/null +++ b/include/mesh/dim3/mass_matrix/mass_matrix.hpp @@ -0,0 +1,92 @@ +#pragma once + +#include "enumerations/dimension.hpp" +#include "mesh/mesh_base.hpp" +#include "specfem_setup.hpp" + +namespace specfem { +namespace mesh { + +/** + * @brief Struct to store mass matrix for a 3D mesh + * + */ +template <> struct mass_matrix { + + constexpr static auto dimension = specfem::dimension::type::dim3; + + int nglob; ///< Number of global nodes + bool is_acoustic; ///< Is acoustic + bool is_elastic; ///< Is elastic + bool is_poroelastic; ///< Is poroelastic + bool has_ocean_load; ///< Has ocean load + + template using View1D = Kokkos::View; + + View1D elastic; ///< Elastic mass matrix + View1D acoustic; ///< Acoustic mass matrix + View1D ocean_load; ///< Ocean load mass matrix + View1D solid_poroelastic; ///< Solid poroelastic mass matrix + View1D fluid_poroelastic; ///< Fluid poroelastic mass matrix + + /** + * @name Constructors + * + */ + ///@{ + /** + * @brief Default constructor initializing an empty struct + * + */ + mass_matrix() = default; + + /** + * @brief Construct a new mass matrix object + * + * @param nglob Number of global nodes + * @param is_acoustic Is acoustic + * @param is_elastic Is elastic + * @param is_poroelastic Is poroelastic + * @param has_ocean_load Has ocean load + * + * @code{.cpp} + * // Example of how to use this constructor + * int nglob = 10; + * bool is_acoustic = true; + * bool is_elastic = false; + * bool is_poroelastic = false; + * bool has_ocean_load = false; + * specfem::mesh::mass_matrix + * mass_matrix(nglob, is_acoustic, is_elastic, is_poroelastic, + * has_ocean_load); + * @endcode + */ + mass_matrix(const int nglob, const bool is_acoustic, const bool is_elastic, + const bool is_poroelastic, const bool has_ocean_load) + : nglob(nglob), is_acoustic(is_acoustic), is_elastic(is_elastic), + is_poroelastic(is_poroelastic), has_ocean_load(has_ocean_load) { + // Only allocate the mass matrix if the simulation includes the respective + // physics. + if (is_acoustic) { + acoustic = + View1D("specfem::mesh::mass_matrix::acoustic", nglob); + } + if (is_elastic) { + elastic = View1D("specfem::mesh::mass_matrix::elastic", nglob); + if (has_ocean_load) { + ocean_load = + View1D("specfem::mesh::mass_matrix::ocean_load", nglob); + } + } + if (is_poroelastic) { + solid_poroelastic = View1D( + "specfem::mesh::mass_matrix::solid_poroelastic", nglob); + fluid_poroelastic = View1D( + "specfem::mesh::mass_matrix::fluid_poroelastic", nglob); + } + }; + ///@} +}; + +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/materials/materials.hpp b/include/mesh/dim3/materials/materials.hpp new file mode 100644 index 000000000..eb59580fe --- /dev/null +++ b/include/mesh/dim3/materials/materials.hpp @@ -0,0 +1,191 @@ +#pragma once +#include "enumerations/dimension.hpp" +#include "mesh/mesh_base.hpp" +#include "specfem_setup.hpp" +#include + +namespace specfem { +namespace mesh { + +/** + * @brief Struct to store materials for a 3D mesh + * + */ +template <> struct materials { + + constexpr static auto dimension = + specfem::dimension::type::dim3; ///< Dimension + ///< type + + int nspec; ///< Number of spectral elements + int ngllx; ///< Number of GLL points in x + int nglly; ///< Number of GLL points in y + int ngllz; ///< Number of GLL points in z + + bool acoustic; ///< Acoustic simulation + bool elastic; ///< Elastic simulation + bool poroelastic; ///< Poroelastic simulation + + template using View4D = Kokkos::View; + template using View5D = Kokkos::View; + + View4D rho; ///< Density rho + View4D kappa; ///< Bulk modulus kappa + View4D mu; ///< Shear modulus mu + + View4D rho_vp; ///< For Stacey ABC + View4D rho_vs; ///< For Stacey ABC + + // Poroelastic properties + View5D poro_rho; ///< Poroelastic density rho + View5D poro_kappa; ///< Poroelastic bulk modulus kappa + View5D poro_perm; ///< Poroelastic permeability perm + + // Poroelastic properties + View4D poro_eta; ///< Poroelastic storage coefficient eta + View4D poro_tort; ///< Poroelastic tortuosity tort + View4D poro_phi; ///< Poroelastic porosity phi + View4D poro_rho_vpI; ///< Poroelastic rho_vpI + View4D poro_rho_vpII; ///< Poroelastic rho_vspII + View4D poro_rho_vsI; ///< Poroelastic rho_vsI + + // Anisotropic properties + bool anisotropic = false; ///< Anisotropic simulation + View4D c11; ///< Anisotropic c11 + View4D c12; ///< Anisotropic c12 + View4D c13; ///< Anisotropic c13 + View4D c14; ///< Anisotropic c14 + View4D c15; ///< Anisotropic c15 + View4D c16; ///< Anisotropic c16 + View4D c22; ///< Anisotropic c22 + View4D c23; ///< Anisotropic c23 + View4D c24; ///< Anisotropic c24 + View4D c25; ///< Anisotropic c25 + View4D c26; ///< Anisotropic c26 + View4D c33; ///< Anisotropic c33 + View4D c34; ///< Anisotropic c34 + View4D c35; ///< Anisotropic c35 + View4D c36; ///< Anisotropic c36 + View4D c44; ///< Anisotropic c44 + View4D c45; ///< Anisotropic c45 + View4D c46; ///< Anisotropic c46 + View4D c55; ///< Anisotropic c55 + View4D c56; ///< Anisotropic c56 + View4D c66; ///< Anisotropic c66 + + /** + * @name Constructors + */ + ///@{ + /** + * @brief Default constructor + * + */ + materials() = default; + + /** + * @brief Construct a new materials object + * + * @param nspec Number of spectral elements + * @param ngllx Number of GLL points in x + * @param nglly Number of GLL points in y + * @param ngllz Number of GLL points in z + * @param acoustic whether the simulation is acoustic + * @param elastic whether the simulation is elastic + * @param poroelastic whether the simulation is poroelastic + * @param anisotropic whether the simulation is anisotropic + */ + materials(const int nspec, const int ngllx, const int nglly, const int ngllz, + const bool acoustic, const bool elastic, const bool poroelastic, + const bool anisotropic) + : nspec(nspec), ngllx(ngllx), nglly(nglly), ngllz(ngllz), + acoustic(acoustic), elastic(elastic), poroelastic(poroelastic), + anisotropic(anisotropic), + rho("specfem::mesh::materials::rho", nspec, ngllx, nglly, ngllz), + kappa("specfem::mesh::materials::kappa", nspec, ngllx, nglly, ngllz), + mu("specfem::mesh::materials::mu", nspec, ngllx, nglly, ngllz) { + if (elastic) { + rho_vp = View4D("specfem::mesh::materials::rho_vp", nspec, + ngllx, nglly, ngllz); + rho_vs = View4D("specfem::mesh::materials::rho_vs", nspec, + ngllx, nglly, ngllz); + + if (anisotropic) { + c11 = View4D("specfem::mesh::materials::c11", nspec, ngllx, + nglly, ngllz); + c12 = View4D("specfem::mesh::materials::c12", nspec, ngllx, + nglly, ngllz); + c13 = View4D("specfem::mesh::materials::c13", nspec, ngllx, + nglly, ngllz); + c14 = View4D("specfem::mesh::materials::c14", nspec, ngllx, + nglly, ngllz); + c15 = View4D("specfem::mesh::materials::c15", nspec, ngllx, + nglly, ngllz); + c16 = View4D("specfem::mesh::materials::c16", nspec, ngllx, + nglly, ngllz); + c22 = View4D("specfem::mesh::materials::c22", nspec, ngllx, + nglly, ngllz); + c23 = View4D("specfem::mesh::materials::c23", nspec, ngllx, + nglly, ngllz); + c24 = View4D("specfem::mesh::materials::c24", nspec, ngllx, + nglly, ngllz); + c25 = View4D("specfem::mesh::materials::c25", nspec, ngllx, + nglly, ngllz); + c26 = View4D("specfem::mesh::materials::c26", nspec, ngllx, + nglly, ngllz); + c33 = View4D("specfem::mesh::materials::c33", nspec, ngllx, + nglly, ngllz); + c34 = View4D("specfem::mesh::materials::c34", nspec, ngllx, + nglly, ngllz); + c35 = View4D("specfem::mesh::materials::c35", nspec, ngllx, + nglly, ngllz); + c36 = View4D("specfem::mesh::materials::c36", nspec, ngllx, + nglly, ngllz); + c44 = View4D("specfem::mesh::materials::c44", nspec, ngllx, + nglly, ngllz); + c45 = View4D("specfem::mesh::materials::c45", nspec, ngllx, + nglly, ngllz); + c46 = View4D("specfem::mesh::materials::c46", nspec, ngllx, + nglly, ngllz); + c55 = View4D("specfem::mesh::materials::c55", nspec, ngllx, + nglly, ngllz); + c56 = View4D("specfem::mesh::materials::c56", nspec, ngllx, + nglly, ngllz); + c66 = View4D("specfem::mesh::materials::c66", nspec, ngllx, + nglly, ngllz); + } + } + if (poroelastic) { + // Hardcoded array sizes + poro_rho = View5D("specfem::mesh::materials::poro_rho", nspec, + 2, ngllx, nglly, ngllz); + poro_kappa = View5D("specfem::mesh::materials::poro_kappa", + nspec, 3, ngllx, nglly, ngllz); + poro_perm = View5D("specfem::mesh::materials::poro_perm", + nspec, 6, ngllx, nglly, ngllz); + poro_eta = View4D("specfem::mesh::materials::poro_eta", nspec, + ngllx, nglly, ngllz); + poro_tort = View4D("specfem::mesh::materials::poro_tort", + nspec, ngllx, nglly, ngllz); + poro_phi = View4D("specfem::mesh::materials::poro_phi", nspec, + ngllx, nglly, ngllz); + poro_rho_vpI = View4D("specfem::mesh::materials::poro_rho_vpI", + nspec, ngllx, nglly, ngllz); + poro_rho_vpII = + View4D("specfem::mesh::materials::poro_rho_vpII", nspec, + ngllx, nglly, ngllz); + poro_rho_vsI = View4D("specfem::mesh::materials::poro_rho_vsI", + nspec, ngllx, nglly, ngllz); + } + }; + + ///@} + + /** + * @brief Print basic information about the materials + * + */ + std::string print(); +}; +} // namespace mesh +} // namespace specfem diff --git a/include/mesh/dim3/mesh.hpp b/include/mesh/dim3/mesh.hpp index f9915afda..2ddbcd5d1 100644 --- a/include/mesh/dim3/mesh.hpp +++ b/include/mesh/dim3/mesh.hpp @@ -1,6 +1,13 @@ #pragma once +#include "boundaries/absorbing_boundary.hpp" +#include "boundaries/boundaries.hpp" +#include "boundaries/free_surface.hpp" #include "coordinates/coordinates.hpp" +#include "coupled_interfaces/coupled_interfaces.hpp" +#include "element_types/element_types.hpp" +#include "mass_matrix/mass_matrix.hpp" +#include "materials/materials.hpp" #include "mesh/dim3/mapping/mapping.hpp" #include "mesh/mesh_base.hpp" #include "parameters/parameters.hpp" @@ -13,6 +20,10 @@ namespace specfem { namespace mesh { +/** + * @brief Struct to store a 3D mesh + * + */ template <> struct mesh { constexpr static auto dimension = @@ -21,60 +32,47 @@ template <> struct mesh { template using View1D = Kokkos::View; // Struct to store all the mesh parameter - specfem::mesh::parameters parameters; + specfem::mesh::parameters parameters; ///< Parameters // Struct to store all the coordinates - specfem::mesh::coordinates coordinates; + specfem::mesh::coordinates coordinates; ///< Coordinates // Struct to store the mapping information - specfem::mesh::mapping mapping; + specfem::mesh::mapping mapping; ///< Mapping // Irregular elements Kokkos - type_real xix_regular, jacobian_regular; - View1D irregular_element_number; + type_real xix_regular, jacobian_regular; ///< Regular xi-xi mapping + View1D irregular_element_number; ///< Irregular elements // Struct to store the partial derivatives - specfem::mesh::partial_derivatives partial_derivatives; - - // - // int npgeo; ///< Total number of spectral element control nodes - // int nspec; ///< Total number of spectral elements - // int nproc; ///< Total number of processors - // specfem::mesh::control_nodes control_nodes; ///< Defines control - // ///< nodes - - // specfem::mesh::parameters parameters; ///< Struct to store - // ///< simulation launch - // ///< parameters (never - // used) - - // specfem::mesh::coupled_interfaces - // coupled_interfaces; ///< Struct to store - // ///< coupled interfaces - - // specfem::mesh::boundaries boundaries; ///< Struct to store - // ///< information at the - // ///< boundaries - - // specfem::mesh::tags tags; ///< Struct to store - // ///< tags for every - // ///< spectral - // ///< element - - // specfem::mesh::elements::tangential_elements - // tangential_nodes; ///< Defines - // ///< tangential - // ///< nodes - // ///< (never - // ///< used) - - // specfem::mesh::elements::axial_elements axial_nodes; ///< - // Defines - // ///< axial - // ///< nodes - // ///< (never - // ///< used) - // specfem::mesh::materials materials; ///< Defines material properties + specfem::mesh::partial_derivatives + partial_derivatives; ///< Partial + ///< derivatives + + // Struct to store element_types + specfem::mesh::element_types elements_types; ///< Element types + + // Mass matrix + specfem::mesh::mass_matrix mass_matrix; ///< Mass matrix + + // Material + specfem::mesh::materials materials; ///< Materials + + // Struct to store the boundaries + specfem::mesh::boundaries boundaries; ///< Boundaries + + // Struct to store the absorbing boundaries + specfem::mesh::absorbing_boundary + absorbing_boundary; ///< Absorbing + ///< boundaries + + // Struct to store the free surface + specfem::mesh::free_surface free_surface; ///< Free surface + + // Struct to store the coupled interfaces + specfem::mesh::coupled_interfaces + coupled_interfaces; ///< Coupled + ///< interfaces /** * @name Constructors @@ -87,30 +85,57 @@ template <> struct mesh { */ mesh(){}; - // mesh(const int npgeo, const int nspec, const int nproc, - // const specfem::mesh::control_nodes - // &control_nodes, - // const specfem::mesh::parameters - // ¶meters, - // const - // specfem::mesh::coupled_interfaces - // &coupled_interfaces, - // const specfem::mesh::boundaries - // &boundaries, - // const specfem::mesh::tags &tags, - // const specfem::mesh::elements::tangential_elements< - // specfem::dimension::type::dim2> &tangential_nodes, - // const specfem::mesh::elements::axial_elements< - // specfem::dimension::type::dim2> &axial_nodes, - // const specfem::mesh::materials &materials) - // : npgeo(npgeo), nspec(nspec), nproc(nproc), - // control_nodes(control_nodes), - // parameters(parameters), coupled_interfaces(coupled_interfaces), - // boundaries(boundaries), tags(tags), - // tangential_nodes(tangential_nodes), axial_nodes(axial_nodes), - // materials(materials){}; - ///@} + /** + * @brief Constructor + * + * @param parameters Struct to store simulation launch parameters + * @param coordinates Struct to store coordinates + * @param mapping Struct to store mapping information + * @param xix_regular Regular xi-xi mapping + * @param jacobian_regular Regular Jacobian + * @param irregular_element_number Kokkos View of irregular elements + * @param partial_derivatives Struct to store partial derivatives + * @param elements_types Struct to store element types + * @param mass_matrix Struct to store mass matrix + * @param materials Struct to store material properties + * @param boundaries Struct to store information at the boundaries + * @param absorbing_boundary Struct to store absorbing boundaries + * @param free_surface Struct to store free surface boundaries + * @param coupled_interfaces Struct to store coupled interfaces + * + * @note This constructor is usually unused, and the mesh is constructed + * using the @ref specfem::IO::read_3d_mesh function. + * + * @see specfem::IO::read_3d_mesh + * + */ + mesh(const specfem::mesh::parameters ¶meters, + const specfem::mesh::coordinates &coordinates, + const specfem::mesh::mapping &mapping, + const type_real xix_regular, const type_real jacobian_regular, + const View1D irregular_element_number, + const specfem::mesh::partial_derivatives &partial_derivatives, + const specfem::mesh::element_types &elements_types, + const specfem::mesh::mass_matrix &mass_matrix, + const specfem::mesh::materials &materials, + const specfem::mesh::boundaries &boundaries, + const specfem::mesh::absorbing_boundary &absorbing_boundary, + const specfem::mesh::free_surface &free_surface, + const specfem::mesh::coupled_interfaces &coupled_interfaces) + : parameters(parameters), coordinates(coordinates), mapping(mapping), + xix_regular(xix_regular), jacobian_regular(jacobian_regular), + irregular_element_number(irregular_element_number), + partial_derivatives(partial_derivatives), + elements_types(elements_types), mass_matrix(mass_matrix), + materials(materials), boundaries(boundaries), + absorbing_boundary(absorbing_boundary), free_surface(free_surface), + coupled_interfaces(coupled_interfaces){}; + + ///@} // Constructors + /** + * @brief Print basic parameters of the mesh if interested + */ std::string print() const; }; } // namespace mesh diff --git a/include/mesh/dim3/parameters/parameters.hpp b/include/mesh/dim3/parameters/parameters.hpp index 6c2addd95..dda6d5658 100644 --- a/include/mesh/dim3/parameters/parameters.hpp +++ b/include/mesh/dim3/parameters/parameters.hpp @@ -65,21 +65,29 @@ template <> struct parameters { int nspec_outer_poroelastic; ///< Number of outer poroelastic SEs // Integer Parameters: coloring - int num_phase_ispec_acoustic; - int num_phase_ispec_elastic; - int num_phase_ispec_poroelastic; - int num_colors_inner_acoustic; - int num_colors_outer_acoustic; - int num_colors_inner_elastic; - int num_colors_outer_elastic; + int num_phase_ispec_acoustic; ///< Number of phase ispec acoustic + int num_phase_ispec_elastic; ///< Number of phase ispec elastic + int num_phase_ispec_poroelastic; ///< Number of phase ispec poroelastic + int num_colors_inner_acoustic; ///< Number of colors inner acoustic + int num_colors_outer_acoustic; ///< Number of colors outer acoustic + int num_colors_inner_elastic; ///< Number of colors inner elastic + int num_colors_outer_elastic; ///< Number of colors outer elastic + /** + * @name Constructors + * + */ + ///@{ /** * @brief Default constructor * */ parameters(){}; - /** Constructor + /** + * + * @brief Construct a new parameters object + * * @param acoustic_simulation Flag for acoustic simulation * @param elastic_simulation Flag for elastic simulation * @param poroelastic_simulation Flag for poroelastic simulation @@ -121,13 +129,13 @@ template <> struct parameters { * @param nspec_inner_poroelastic Number of inner poroelastic SEs * @param nspec_outer_poroelastic Number of outer poroelastic SEs * - * @param num_phase_ispec_acoustic - * @param num_phase_ispec_elastic - * @param num_phase_ispec_poroelastic - * @param num_colors_inner_acoustic - * @param num_colors_outer_acoustic - * @param num_colors_inner_elastic - * @param num_colors_outer_elastic + * @param num_phase_ispec_acoustic Number of phase ispec acoustic + * @param num_phase_ispec_elastic Number of phase ispec elastic + * @param num_phase_ispec_poroelastic Number of phase ispec poroelastic + * @param num_colors_inner_acoustic Number of colors inner acoustic + * @param num_colors_outer_acoustic Number of colors outer acoustic + * @param num_colors_inner_elastic Number of colors inner elastic + * @param num_colors_outer_elastic Number of colors outer elastic * */ parameters( @@ -187,7 +195,13 @@ template <> struct parameters { num_colors_inner_elastic(num_colors_inner_elastic), num_colors_outer_elastic(num_colors_outer_elastic){}; - void print() const; + ///@} + + /** + * @brief Print basic information about the parameters + * + */ + std::string print() const; }; } // namespace mesh diff --git a/include/mesh/dim3/partial_derivatives/partial_derivatives.hpp b/include/mesh/dim3/partial_derivatives/partial_derivatives.hpp index 05ce270cd..b803d0336 100644 --- a/include/mesh/dim3/partial_derivatives/partial_derivatives.hpp +++ b/include/mesh/dim3/partial_derivatives/partial_derivatives.hpp @@ -1,39 +1,56 @@ #pragma once #include "enumerations/dimension.hpp" +#include "mesh/mesh_base.hpp" #include "specfem_setup.hpp" #include namespace specfem { namespace mesh { -template struct partial_derivatives; - +/** + * @brief Struct to store partial derivatives for a 3D mesh + * + */ template <> struct partial_derivatives { constexpr static auto dimension = specfem::dimension::type::dim3; using LocalView = Kokkos::View; // Parameters - int nspec; - int ngllx; - int nglly; - int ngllz; - - LocalView xix; - LocalView xiy; - LocalView xiz; - LocalView etax; - LocalView etay; - LocalView etaz; - LocalView gammax; - LocalView gammay; - LocalView gammaz; - LocalView jacobian; - - // Constructors + int nspec; ///< Number of spectral elements + int ngllx; ///< Number of GLL points in x + int nglly; ///< Number of GLL points in y + int ngllz; ///< Number of GLL points in z + + LocalView xix; ///< 4D Kokkos::view of type real for xix + LocalView xiy; ///< 4D Kokkos::view of type real for xiy + LocalView xiz; ///< 4D Kokkos::view of type real for xiz + LocalView etax; ///< 4D Kokkos::view of type real for etax + LocalView etay; ///< 4D Kokkos::view of type real for etay + LocalView etaz; ///< 4D Kokkos::view of type real for etaz + LocalView gammax; ///< 4D Kokkos::view of type real for gammax + LocalView gammay; ///< 4D Kokkos::view of type real for gammay + LocalView gammaz; ///< 4D Kokkos::view of type real for gammaz + LocalView jacobian; ///< 4D Kokkos::view of type real for jacobian + + /** + * @name Constructors + */ + ///@{ + /** + * @brief Default constructor + * + */ partial_derivatives(){}; // Default constructor - // Constructor to initialize the coordinates + /** + * @brief Construct a new partial derivatives object + * + * @param nspec Number of spectral elements + * @param ngllx Number of GLL points in x + * @param nglly Number of GLL points in y + * @param ngllz Number of GLL points in z + */ partial_derivatives(int nspec, int ngllx, int nglly, int ngllz) : nspec(nspec), ngllx(ngllx), nglly(nglly), ngllz(ngllz), xix("xix", nspec, ngllx, nglly, ngllz), @@ -47,9 +64,24 @@ template <> struct partial_derivatives { gammaz("gammaz", nspec, ngllx, nglly, ngllz), jacobian("jacobian", nspec, ngllx, nglly, ngllz){}; - void print() const; + ///@} + + /** + * @brief Print basic information on the partial derivatives + * + */ + std::string print() const; - void print(int ispec, int igllx, int iglly, int igllz) const; + /** + * @brief Print the partial derivatives for a given spectral element and + * GLL point + * + * @param ispec Spectral element index + * @param igllx GLL point index in x + * @param iglly GLL point index in y + * @param igllz GLL point index in z + */ + std::string print(int ispec, int igllx, int iglly, int igllz) const; }; } // namespace mesh diff --git a/include/mesh/mesh_base.hpp b/include/mesh/mesh_base.hpp index 2d9c2af7a..b2ea63c91 100644 --- a/include/mesh/mesh_base.hpp +++ b/include/mesh/mesh_base.hpp @@ -8,57 +8,116 @@ namespace mesh { /** * @brief Struct to store information about the mesh read from the database * + * @tparam DimensionType Dimension type */ template struct mesh; /** * @brief Struct to store general parameters for the mesh * + * @tparam DimensionType Dimension type */ template struct parameters; +/** + * @brief Struct to store materials + * + * @tparam DimensionType Dimension type + */ +template struct materials; + +/** + * @brief Struct to store coordinates + * + * @tparam DimensionType Dimension type + */ +template struct coordinates; + +/** + * @brief Struct to store mapping + * + * @tparam DimensionType Dimension type + */ +template struct mapping; + +/** + * @brief Struct to store partial derivatives + * + * @tparam DimensionType Dimension type + */ +template struct partial_derivatives; + /** * @brief Struct to store control nodes * + * @tparam DimensionType Dimension type */ template struct control_nodes; /** * @brief Struct to store coupled interfaces * + * @tparam DimensionType Dimension type */ template struct coupled_interfaces; /** * @brief Struct to store boundaries * + * @tparam DimensionType Dimension type */ template struct boundaries; /** * @brief Struct to store absorbing boundaries * + * @tparam DimensionType Dimension type */ template struct absorbing_boundary; /** * @brief Struct to store acoustic free surface boundaries * + * @tparam DimensionType Dimension type */ template struct acoustic_free_surface; +/** + * @brief Struct to store acoustic free surface boundaries + * + * @tparam DimensionType Dimension type + */ +template struct free_surface; + /** * @brief Struct to store forcing boundaries * + * @tparam DimensionType Dimension type */ template struct forcing_boundary; /** * @brief Struct to store tags for every spectral element * + * @tparam DimensionType Dimension type */ template struct tags; +/** + * @brief Struct to store mass_matrices + * + * @tparam DimensionType Dimension type + */ +template struct mass_matrix; + +/** + * @brief Struct to store element_types (Acoustic, lastic, Poroelastic) + * + * @tparam DimensionType Dimension type + * + */ +template struct element_types; + namespace elements { template struct axial_elements; diff --git a/pyproject.toml b/pyproject.toml index 79ac4e363..cf0401013 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,7 +29,8 @@ docs = [ "sphinx-copybutton~=0.4.0", "furo~=2024.8.6", "breathe~=4.35.0", - "sphinx-sitemap~=2.2.0" + "sphinx-sitemap~=2.2.0", + "sphinx_design" ] [tool.scikit-build] diff --git a/src/IO/fortranio/fortran_io.cpp b/src/IO/fortranio/fortran_io.cpp index a86b470a6..0b76e3260 100644 --- a/src/IO/fortranio/fortran_io.cpp +++ b/src/IO/fortranio/fortran_io.cpp @@ -1,4 +1,5 @@ #include "IO/fortranio/fortran_io.hpp" +#include "IO/fortranio/fortran_io.tpp" #include #include #include @@ -90,3 +91,19 @@ void specfem::IO::fortran_read_value(std::string *value, std::ifstream &stream, boost::algorithm::trim(*value); return; } + +template <> +void specfem::IO::fortran_read_value(std::vector *value, + std::ifstream &stream, + int &buffer_length) { + int nsize = value->size(); + std::vector &rvalue = *value; + for (int i = 0; i < nsize; i++) { + // Create a temporary bool variable to hold the value + bool temp_bool; + specfem::IO::fortran_read_value(&temp_bool, stream, buffer_length); + // Assign the temporary value to the vector element + rvalue[i] = temp_bool; + } + return; +} diff --git a/src/IO/mesh/impl/fortran/dim2/read_material_properties.cpp b/src/IO/mesh/impl/fortran/dim2/read_material_properties.cpp index b40cfa317..f09aa3f67 100644 --- a/src/IO/mesh/impl/fortran/dim2/read_material_properties.cpp +++ b/src/IO/mesh/impl/fortran/dim2/read_material_properties.cpp @@ -1,5 +1,6 @@ #include "IO/mesh/impl/fortran/dim2/read_material_properties.hpp" #include "IO/fortranio/interface.hpp" +#include "enumerations/dimension.hpp" #include "enumerations/interface.hpp" #include "mesh/mesh.hpp" #include "specfem_mpi/interface.hpp" @@ -20,20 +21,25 @@ struct input_holder { int n, indic; }; -std::vector read_materials( +std::vector::material_specification> +read_materials( std::ifstream &stream, const int numat, - specfem::mesh::materials::material &elastic_isotropic, - specfem::mesh::materials::material &acoustic_isotropic, - specfem::mesh::materials::material - &elastic_anisotropic, + specfem::mesh::materials::material< + elastic, isotropic> &elastic_isotropic, + specfem::mesh::materials::material< + acoustic, isotropic> &acoustic_isotropic, + specfem::mesh::materials::material< + elastic, anisotropic> &elastic_anisotropic, const specfem::MPI::MPI *mpi) { input_holder read_values; std::ostringstream message; - std::vector index_mapping( - numat); + std::vector::material_specification> + index_mapping(numat); message << "Material systems:\n" << "------------------------------"; @@ -104,10 +110,11 @@ std::vector read_materials( l_acoustic_isotropic.push_back(acoustic_isotropic_holder); - index_mapping[i] = specfem::mesh::materials::material_specification( - specfem::element::medium_tag::acoustic, - specfem::element::property_tag::isotropic, index_acoustic_isotropic, - read_values.n - 1); + index_mapping[i] = specfem::mesh:: + materials::material_specification( + specfem::element::medium_tag::acoustic, + specfem::element::property_tag::isotropic, + index_acoustic_isotropic, read_values.n - 1); index_acoustic_isotropic++; @@ -128,10 +135,11 @@ std::vector read_materials( l_elastic_isotropic.push_back(elastic_isotropic_holder); - index_mapping[i] = specfem::mesh::materials::material_specification( - specfem::element::medium_tag::elastic, - specfem::element::property_tag::isotropic, index_elastic_isotropic, - read_values.n - 1); + index_mapping[i] = specfem::mesh:: + materials::material_specification( + specfem::element::medium_tag::elastic, + specfem::element::property_tag::isotropic, + index_elastic_isotropic, read_values.n - 1); index_elastic_isotropic++; } @@ -159,10 +167,11 @@ std::vector read_materials( l_elastic_anisotropic.push_back(elastic_anisotropic_holder); - index_mapping[i] = specfem::mesh::materials::material_specification( - specfem::element::medium_tag::elastic, - specfem::element::property_tag::anisotropic, - index_elastic_anisotropic, read_values.n - 1); + index_mapping[i] = specfem::mesh:: + materials::material_specification( + specfem::element::medium_tag::elastic, + specfem::element::property_tag::anisotropic, + index_elastic_anisotropic, read_values.n - 1); index_elastic_anisotropic++; @@ -175,25 +184,29 @@ std::vector read_materials( l_elastic_anisotropic.size() == numat); - elastic_isotropic = specfem::mesh::materials::material( - l_elastic_isotropic.size(), l_elastic_isotropic); + elastic_isotropic = + specfem::mesh::materials::material< + elastic, isotropic>(l_elastic_isotropic.size(), l_elastic_isotropic); elastic_anisotropic = - specfem::mesh::materials::material( - l_elastic_anisotropic.size(), l_elastic_anisotropic); + specfem::mesh::materials::material< + elastic, anisotropic>(l_elastic_anisotropic.size(), + l_elastic_anisotropic); - acoustic_isotropic = specfem::mesh::materials::material( - l_acoustic_isotropic.size(), l_acoustic_isotropic); + acoustic_isotropic = + specfem::mesh::materials::material< + acoustic, isotropic>(l_acoustic_isotropic.size(), + l_acoustic_isotropic); return index_mapping; } void read_material_indices( std::ifstream &stream, const int nspec, const int numat, - const std::vector - &index_mapping, - const specfem::kokkos::HostView1d< - specfem::mesh::materials::material_specification> + const std::vector::material_specification> &index_mapping, + const specfem::kokkos::HostView1d::material_specification> material_index_mapping, const specfem::kokkos::HostView2d knods, const specfem::MPI::MPI *mpi) { @@ -230,14 +243,15 @@ void read_material_indices( return; } -specfem::mesh::materials +specfem::mesh::materials specfem::IO::mesh::impl::fortran::dim2::read_material_properties( std::ifstream &stream, const int numat, const int nspec, const specfem::kokkos::HostView2d knods, const specfem::MPI::MPI *mpi) { // Create materials instances - specfem::mesh::materials materials(nspec, numat); + specfem::mesh::materials materials(nspec, + numat); // Read material properties auto index_mapping = read_materials( diff --git a/src/IO/mesh/impl/fortran/dim3/mesh.cpp b/src/IO/mesh/impl/fortran/dim3/mesh.cpp index be961ddb1..12d4d7c0a 100644 --- a/src/IO/mesh/impl/fortran/dim3/mesh.cpp +++ b/src/IO/mesh/impl/fortran/dim3/mesh.cpp @@ -2,12 +2,6 @@ #include "mesh/mesh.hpp" #include "IO/fortranio/interface.hpp" #include "IO/interface.hpp" -#include "IO/mesh/impl/fortran/dim2/read_boundaries.hpp" -#include "IO/mesh/impl/fortran/dim2/read_elements.hpp" -#include "IO/mesh/impl/fortran/dim2/read_interfaces.hpp" -#include "IO/mesh/impl/fortran/dim2/read_material_properties.hpp" -#include "IO/mesh/impl/fortran/dim2/read_mesh_database.hpp" -#include "IO/mesh/impl/fortran/dim2/read_parameters.hpp" #include "IO/mesh/impl/fortran/dim3/interface.hpp" #include "enumerations/interface.hpp" #include "kokkos_abstractions.h" @@ -18,16 +12,47 @@ // External/Standard Libraries #include #include +#include +#include #include #include +#include +#include #include #include +template +std::string try_print_medium_element( + const specfem::mesh::element_types + &elements_types, + int index) { + + std::ostringstream message; + try { + message << elements_types.print(index); + } catch (std::runtime_error &e) { + message << e.what(); + } catch (...) { + message << "Unknown exception caught in try_print_medium_element" + << ".\n"; + } + return message.str(); +} + specfem::mesh::mesh specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, const std::string mesh_databases_file, const specfem::MPI::MPI *mpi) { + // Creating aliases for Checking functions + using specfem::IO::mesh::impl::fortran::dim3::check_read_test_value; + using specfem::IO::mesh::impl::fortran::dim3::check_values; + + // Creating aliases for Array Reading functions + using specfem::IO::mesh::impl::fortran::dim3::try_read_array; + using specfem::IO::mesh::impl::fortran::dim3::try_read_index_array; + using specfem::IO::mesh::impl::fortran::dim3::try_read_line; + // Declaring empty mesh objects specfem::mesh::mesh mesh; @@ -47,11 +72,14 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, throw; } - mesh.parameters.print(); - // Mesh Parameters are populated, closing the parameters file. stream.close(); +#ifndef NDEBUG + // Print the parameters + mpi->cout(mesh.parameters.print()); +#endif + // Open the database file stream.open(mesh_databases_file); @@ -61,49 +89,15 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, int nspec, nglob, nspec_irregular; - try { - specfem::IO::fortran_read_line(stream, &nspec); - specfem::IO::fortran_read_line(stream, &nglob); - specfem::IO::fortran_read_line(stream, &nspec_irregular); - - std::tie(nspec, nglob, nspec_irregular) = - std::make_tuple(nspec, nglob, nspec_irregular); - - } catch (const std::exception &e) { - std::ostringstream error_message; - error_message - << "Error reading nspec, nglob, nspec_irregular from database file: " - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - } + try_read_line("read_nspec", stream, &nspec); + try_read_line("read_nglob", stream, &nglob); + try_read_line("read_nspec_irregular", stream, &nspec_irregular); - if (nspec != mesh.parameters.nspec) { - std::ostringstream error_message; - error_message << "Database nspec not equal to mesh parameters nspec. " - << "Database nspec = " << nspec - << " Mesh parameters nspec = " << mesh.parameters.nspec << "(" - << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } - - if (nglob != mesh.parameters.nglob) { - std::ostringstream error_message; - error_message << "Database nglob not equal to mesh parameters nglob. " - << "Database nglob = " << nglob - << " Mesh parameters nglob = " << mesh.parameters.nglob << "(" - << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } - - if (nspec_irregular != mesh.parameters.nspec_irregular) { - std::ostringstream error_message; - error_message << "Database nspec_irregular not equal to mesh parameters " - "nspec_irregular. " - << "Database nspec_irregular = " << nspec_irregular - << " Mesh parameters nspec_irregular = " - << mesh.parameters.nspec_irregular << "(" << __FILE__ << ":" - << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } + // Check values + check_values("nspec", nspec, mesh.parameters.nspec); + check_values("nglob", nglob, mesh.parameters.nglob); + check_values("nspec_irregular", nspec_irregular, + mesh.parameters.nspec_irregular); // Create the mapping object mesh.mapping = specfem::mesh::mapping( @@ -112,20 +106,14 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, mesh.parameters.nglly, mesh.parameters.ngllz); // Reading the mapping from the database file. - try { - specfem::IO::mesh::impl::fortran::dim3::read_index_array( - stream, mesh.mapping.ibool); - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading ibool from database file: " << e.what() - << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } + try_read_index_array("read_ibool", stream, mesh.mapping.ibool); +#ifndef NDEBUG // Print Mapping parameters and the first spectral element - mesh.mapping.print(); - mesh.mapping.print(0); - mesh.mapping.print(mesh.parameters.nspec - 1); + mpi->cout(mesh.mapping.print()); + mpi->cout(mesh.mapping.print(0)); + mpi->cout(mesh.mapping.print(mesh.parameters.nspec - 1)); +#endif // Create the coordinates object mesh.coordinates = specfem::mesh::coordinates( @@ -133,68 +121,43 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, mesh.parameters.nglly, mesh.parameters.ngllz); // Reading the coordinates from the database file. - try { - specfem::IO::mesh::impl::fortran::dim3::read_xyz(stream, mesh.coordinates, - mpi); - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading xyz from database file: " << e.what() << "(" - << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } + try_read_array("read_x", stream, mesh.coordinates.x); + try_read_array("read_y", stream, mesh.coordinates.y); + try_read_array("read_z", stream, mesh.coordinates.z); +#ifndef NDEBUG // Print Coordinates parameters and the first global node - mesh.coordinates.print(); - mesh.coordinates.print(0); - mesh.coordinates.print(mesh.parameters.nglob - 1); + mpi->cout(mesh.coordinates.print()); + mpi->cout(mesh.coordinates.print(0)); + mpi->cout(mesh.coordinates.print(mesh.parameters.nglob - 1)); +#endif // Initialize the partial derivatives object mesh.irregular_element_number = decltype(mesh.irregular_element_number)( "irregular_element_number", mesh.parameters.nspec_irregular); // Read Irregular elements - try { - specfem::IO::mesh::impl::fortran::dim3::read_index_array( - stream, mesh.irregular_element_number); - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message - << "Error reading irregular_element_number from database file: " - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } + try_read_index_array("read_irregular_element_number", stream, + mesh.irregular_element_number); // Read the partial derivatives (only two CUSTOM_REALs) - try { - fortran_read_line(stream, &mesh.xix_regular); - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading xix_regular from database file: " - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } - - // Read the partial derivatives (only two CUSTOM_REALs) - try { - fortran_read_line(stream, &mesh.jacobian_regular); - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading jacobian_regular from database file: " - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); - } + try_read_line("read_xix_regular", stream, &mesh.xix_regular); + try_read_line("read_jacobian_regular", stream, &mesh.jacobian_regular); +#ifndef NDEBUG // Print the first and last irregular element - std::cout << "First irregular element: " << mesh.irregular_element_number(0) - << std::endl; - std::cout << "Last irregular element: " - << mesh.irregular_element_number(mesh.parameters.nspec_irregular - - 1) - << std::endl; + std::ostringstream message; + message << "First irregular element: " << mesh.irregular_element_number(0) + << "\n"; + message << "Last irregular element: " + << mesh.irregular_element_number(mesh.parameters.nspec_irregular - 1) + << "\n"; // Print xix and jacobian - std::cout << "xix_regular: " << mesh.xix_regular << std::endl; - std::cout << "jacobian_regular: " << mesh.jacobian_regular << std::endl; + message << "xix_regular: " << mesh.xix_regular << "\n"; + message << "jacobian_regular: " << mesh.jacobian_regular << "\n"; + mpi->cout(message.str()); +#endif // Create the partial derivatives object mesh.partial_derivatives = @@ -203,19 +166,375 @@ specfem::IO::read_3d_mesh(const std::string mesh_parameters_file, mesh.parameters.ngllz); // Reading the partial derivatives from the database file. - try { - specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( - stream, mesh.partial_derivatives, mpi); - } catch (std::runtime_error &e) { - std::ostringstream error_message; - error_message << "Error reading partial derivatives from database file: " - << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; - throw std::runtime_error(error_message.str()); + try_read_array("read_xi_x", stream, mesh.partial_derivatives.xix); + try_read_array("read_xi_y", stream, mesh.partial_derivatives.xiy); + try_read_array("read_xi_z", stream, mesh.partial_derivatives.xiz); + try_read_array("read_eta_x", stream, mesh.partial_derivatives.etax); + try_read_array("read_eta_y", stream, mesh.partial_derivatives.etay); + try_read_array("read_eta_z", stream, mesh.partial_derivatives.etaz); + try_read_array("read_gamma_x", stream, mesh.partial_derivatives.gammax); + try_read_array("read_gamma_y", stream, mesh.partial_derivatives.gammay); + try_read_array("read_gamma_z", stream, mesh.partial_derivatives.gammaz); + try_read_array("read_jacobian", stream, mesh.partial_derivatives.jacobian); + +#ifndef NDEBUG + // Print Partial Derivatives parameters and the first spectral element + mpi->cout(mesh.partial_derivatives.print()); + mpi->cout(mesh.partial_derivatives.print(0, 0, 0, 0)); +#endif + + // Marker that should be 10000 + check_read_test_value(stream, 10000); + + // Create material object + mesh.materials = specfem::mesh::materials( + mesh.parameters.nspec, mesh.parameters.ngllx, mesh.parameters.nglly, + mesh.parameters.ngllz, mesh.parameters.acoustic_simulation, + mesh.parameters.elastic_simulation, + mesh.parameters.poroelastic_simulation, mesh.parameters.anisotropy); + + try_read_array("read_kappa", stream, mesh.materials.kappa); + try_read_array("read_mu", stream, mesh.materials.mu); + +#ifndef NDEBUG + // Print the materials + mpi->cout(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, nacoustic, nelastic, nporoelastic); + + // 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(ispec_type); + +#ifndef NDEBUG + // Print the element types + mpi->cout(mesh.elements_types.print()); + mpi->cout(mesh.elements_types.print(0)); + + // Print elements first element of each category + 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 + check_read_test_value(stream, 9999); + + // Intialize the mass matrix object + mesh.mass_matrix = specfem::mesh::mass_matrix( + mesh.parameters.nglob, mesh.parameters.acoustic_simulation, + mesh.parameters.elastic_simulation, + mesh.parameters.poroelastic_simulation, + mesh.parameters.approximate_ocean_load); + + // Read the acoustic mass matrix if acoustic simulation + if (mesh.parameters.acoustic_simulation) { + try_read_array("read_acoustic_mass_matrix", stream, + mesh.mass_matrix.acoustic); } - // Print Partial Derivatives parameters and the first spectral element - mesh.partial_derivatives.print(); - mesh.partial_derivatives.print(0, 0, 0, 0); + // Read the density rho + try_read_array("read_rho", stream, mesh.materials.rho); + + // Read test value 9998 + check_read_test_value(stream, 9998); + + // Read the elastic mass matrix if elastic simulation + if (mesh.parameters.elastic_simulation) { + try_read_array("read_elastic_mass_matrix", stream, + mesh.mass_matrix.elastic); + + // Read the stacey boundary values + try_read_array("read_stacey_boundary_values", stream, + mesh.materials.rho_vp); + try_read_array("read_stacey_boundary_values", stream, + mesh.materials.rho_vs); + } + + // Read test value 9997 + check_read_test_value(stream, 9997); + + // If simulation poroelastic + if (mesh.parameters.poroelastic_simulation) { + try_read_array("read_poroelastic_mass_matrix_solid", stream, + mesh.mass_matrix.solid_poroelastic); + try_read_array("read_poroelastic_mass_matrix_fluid", stream, + mesh.mass_matrix.fluid_poroelastic); + + // read the poroelastic material properties + try_read_array("read_poroelastic_rho", stream, mesh.materials.poro_rho); + try_read_array("read_poroelastic_kappa", stream, mesh.materials.poro_kappa); + try_read_array("read_poroelastic_eta", stream, mesh.materials.poro_eta); + try_read_array("read_poroelastic_tort", stream, mesh.materials.poro_tort); + try_read_array("read_poroelastic_perm", stream, mesh.materials.poro_perm); + try_read_array("read_poroelastic_phi", stream, mesh.materials.poro_phi); + try_read_array("read_poroelastic_rho_vpI", stream, + mesh.materials.poro_rho_vpI); + try_read_array("read_poroelastic_rho_vpII", stream, + mesh.materials.poro_rho_vpII); + try_read_array("read_poroelastic_rho_vsI", stream, + mesh.materials.poro_rho_vsI); + } + + // Read test value 9996 + check_read_test_value(stream, 9996); + + // Read number of absorbing boundaries + int num_abs_boundary_faces; + try_read_line("num_abs_boundary_faces", stream, &num_abs_boundary_faces); + + // Check whether the number of boundaries faces is equal to the number of + // faces in the mesh parameters + check_values("num_abs_boundary_faces", num_abs_boundary_faces, + mesh.parameters.num_abs_boundary_faces); + + // if there are absorbing boundaries create the absorbing boundary object + // and read the absorbing boundaries + if (num_abs_boundary_faces > 0) { + mesh.absorbing_boundary = + specfem::mesh::absorbing_boundary( + mesh.parameters.nglob, mesh.parameters.num_abs_boundary_faces, + mesh.parameters.ngllsquare, mesh.parameters.acoustic_simulation, + mesh.parameters.elastic_simulation, mesh.parameters.nspec2D_xmin, + mesh.parameters.nspec2D_xmax, mesh.parameters.nspec2D_ymin, + mesh.parameters.nspec2D_ymax, mesh.parameters.nspec2D_bottom, + mesh.parameters.nspec2D_top); + + // Read the absorbing ispec + try_read_index_array("abs_boundary_ispec", stream, + mesh.absorbing_boundary.ispec); + try_read_index_array("abs_bounary_ijk", stream, + mesh.absorbing_boundary.ijk); + try_read_array("abs_boundary_jacobian2Dw", stream, + mesh.absorbing_boundary.jacobian2Dw); + try_read_array("abs_boundary_normal", stream, + mesh.absorbing_boundary.normal); + + // Read the absorbing mass matrix elastic + if (mesh.parameters.elastic_simulation) { + try_read_array("abs_boundary_mass_elastic.mass", stream, + mesh.absorbing_boundary.mass_elastic.x); + try_read_array("abs_boundary_mass_elastic.mass", stream, + mesh.absorbing_boundary.mass_elastic.y); + try_read_array("abs_boundary_mass_elastic.mass", stream, + mesh.absorbing_boundary.mass_elastic.z); + } + + // Read the absorbing mass matrix acoustic + if (mesh.parameters.acoustic_simulation) { + try_read_array("abs_boundary_mass_acoustic.mass", stream, + mesh.absorbing_boundary.mass_acoustic.mass); + } +#ifndef NDEBUG + // Print the absorbing boundaries + mpi->cout(mesh.absorbing_boundary.print()); +#endif + } + + // Read test value 9995 + check_read_test_value(stream, 9995); + + // Read the number of boundaries + int nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, + nspec2D_top; + + try_read_line("nspec2D_xmin", stream, &nspec2D_xmin); + try_read_line("nspec2D_xmax", stream, &nspec2D_xmax); + try_read_line("nspec2D_ymin", stream, &nspec2D_ymin); + try_read_line("nspec2D_ymax", stream, &nspec2D_ymax); + try_read_line("nspec2D_bottom", stream, &nspec2D_bottom); + try_read_line("nspec2D_top", stream, &nspec2D_top); + +#ifndef NDEBUG + // Print the number of boundaries + message = std::ostringstream(); + message << "nspec2D_xmin: " << nspec2D_xmin << "\n"; + message << "nspec2D_xmax: " << nspec2D_xmax << "\n"; + message << "nspec2D_ymin: " << nspec2D_ymin << "\n"; + message << "nspec2D_ymax: " << nspec2D_ymax << "\n"; + message << "nspec2D_bottom: " << nspec2D_bottom << "\n"; + message << "nspec2D_top: " << nspec2D_top << "\n"; + mpi->cout(message.str()); +#endif + + // Check values + check_values("nspec2D_xmin", nspec2D_xmin, mesh.parameters.nspec2D_xmin); + check_values("nspec2D_xmax", nspec2D_xmax, mesh.parameters.nspec2D_xmax); + check_values("nspec2D_ymin", nspec2D_ymin, mesh.parameters.nspec2D_ymin); + check_values("nspec2D_ymax", nspec2D_ymax, mesh.parameters.nspec2D_ymax); + check_values("nspec2D_bottom", nspec2D_bottom, + mesh.parameters.nspec2D_bottom); + check_values("nspec2D_top", nspec2D_top, mesh.parameters.nspec2D_top); + + if (mesh.parameters.num_abs_boundary_faces > 0) { + + if (nspec2D_xmin > 0) { + try_read_index_array("ibelm_xmin", stream, + mesh.absorbing_boundary.ibelm_xmin); + } + if (nspec2D_xmax > 0) { + try_read_index_array("ibelm_xmax", stream, + mesh.absorbing_boundary.ibelm_xmax); + } + if (nspec2D_ymin > 0) { + try_read_index_array("ibelm_ymin", stream, + mesh.absorbing_boundary.ibelm_ymin); + } + if (nspec2D_ymax > 0) { + try_read_index_array("ibelm_ymax", stream, + mesh.absorbing_boundary.ibelm_ymax); + } + if (nspec2D_bottom > 0) { + try_read_index_array("ibelm_bottom", stream, + mesh.absorbing_boundary.ibelm_bottom); + } + if (nspec2D_top > 0) { + try_read_index_array("ibelm_top", stream, + mesh.absorbing_boundary.ibelm_top); + } + // Print the abosrbing boundaries +#ifndef NDEBUG + mpi->cout(mesh.absorbing_boundary.print()); +#endif + } + + // Read the number of free surface faces + int num_free_surface_faces; + try_read_line("num_free_surface_faces", stream, &num_free_surface_faces); + check_values("num_free_surface_faces", num_free_surface_faces, + mesh.parameters.num_free_surface_faces); + + // Create the free surface object + mesh.free_surface = + specfem::mesh::free_surface( + mesh.parameters.num_free_surface_faces, mesh.parameters.ngllsquare); + + // Read the free surface ispec + try_read_index_array("free_surface_ispec", stream, mesh.free_surface.ispec); + try_read_index_array("free_surface_ijk", stream, mesh.free_surface.ijk); + try_read_array("free_surface_jacobian2Dw", stream, + mesh.free_surface.jacobian2Dw); + try_read_array("free_surface_normal", stream, mesh.free_surface.normal); + +#ifndef NDEBUG + // Print the free surface + mpi->cout(mesh.free_surface.print()); +#endif + + // Create the coupled interfaces object + mesh.coupled_interfaces = + specfem::mesh::coupled_interfaces( + mesh.parameters.num_coupling_ac_el_faces, + mesh.parameters.num_coupling_ac_po_faces, + mesh.parameters.num_coupling_el_po_faces, mesh.parameters.ngllsquare); + + // Read the num_coupling_ac_el_faces + int num_coupling_ac_el_faces; + try_read_line("num_coupling_ac_el_faces", stream, &num_coupling_ac_el_faces); + check_values("num_coupling_ac_el_faces", num_coupling_ac_el_faces, + mesh.parameters.num_coupling_ac_el_faces); + + // Read the the coupling if the number is greater than 0 + if (mesh.coupled_interfaces.acoustic_elastic) { + try_read_index_array( + "coupling_ac_el_ispec", stream, + mesh.coupled_interfaces.acoustic_elastic_interface.ispec); + try_read_index_array( + "coupling_ac_el_ijk", stream, + mesh.coupled_interfaces.acoustic_elastic_interface.ijk); + try_read_array( + "coupling_ac_el_jacobian2Dw", stream, + mesh.coupled_interfaces.acoustic_elastic_interface.jacobian2Dw); + try_read_array("coupling_ac_el_normal", stream, + mesh.coupled_interfaces.acoustic_elastic_interface.normal); + } + + // Read the num_coupling_ac_po_faces + int num_coupling_ac_po_faces; + try_read_line("num_coupling_ac_po_faces", stream, &num_coupling_ac_po_faces); + + // Check the values + check_values("num_coupling_ac_po_faces", num_coupling_ac_po_faces, + mesh.parameters.num_coupling_ac_po_faces); + + // Read the the coupling if the number is greater than 0 + if (mesh.coupled_interfaces.acoustic_poroelastic) { + try_read_index_array( + "coupling_ac_po_ispec", stream, + mesh.coupled_interfaces.acoustic_poroelastic_interface.ispec); + try_read_index_array( + "coupling_ac_po_ijk", stream, + mesh.coupled_interfaces.acoustic_poroelastic_interface.ijk); + try_read_array( + "coupling_ac_po_jacobian2Dw", stream, + mesh.coupled_interfaces.acoustic_poroelastic_interface.jacobian2Dw); + try_read_array( + "coupling_ac_po_normal", stream, + mesh.coupled_interfaces.acoustic_poroelastic_interface.normal); + } + + // Read the num_coupling_el_po_faces + int num_coupling_el_po_faces; + try_read_line("num_coupling_el_po_faces", stream, &num_coupling_el_po_faces); + + // Check the values + check_values("num_coupling_el_po_faces", num_coupling_el_po_faces, + mesh.parameters.num_coupling_el_po_faces); + + // Read the the coupling if the number is greater than 0 + if (mesh.coupled_interfaces.elastic_poroelastic) { + try_read_index_array( + "coupling_el_po_ispec", stream, + mesh.coupled_interfaces.elastic_poroelastic_interface.ispec); + try_read_index_array( + "coupling_el_po_ijk", stream, + mesh.coupled_interfaces.elastic_poroelastic_interface.ijk); + try_read_array( + "coupling_el_po_jacobian2Dw", stream, + mesh.coupled_interfaces.elastic_poroelastic_interface.jacobian2Dw); + try_read_array( + "coupling_el_po_normal", stream, + mesh.coupled_interfaces.elastic_poroelastic_interface.normal); + try_read_index_array( + "coupling_po_el_ispec", stream, + mesh.coupled_interfaces.poroelastic_elastic_interface.ispec); + try_read_index_array( + "coupling_po_el_ijk", stream, + mesh.coupled_interfaces.poroelastic_elastic_interface.ijk); + try_read_array( + "coupling_po_el_jacobian2Dw", stream, + mesh.coupled_interfaces.poroelastic_elastic_interface.jacobian2Dw); + try_read_array( + "coupling_po_el_normal", stream, + mesh.coupled_interfaces.poroelastic_elastic_interface.normal); + } + +#ifndef NDEBUG + // Print the interfaces + mpi->cout(mesh.coupled_interfaces.print()); +#endif + + // Read test value 9997 + check_read_test_value(stream, 9997); + + // Final print with basic information + mpi->cout(mesh.print()); stream.close(); diff --git a/src/IO/mesh/impl/fortran/dim3/read_parameters.cpp b/src/IO/mesh/impl/fortran/dim3/read_parameters.cpp index 5599722b9..909febaf1 100644 --- a/src/IO/mesh/impl/fortran/dim3/read_parameters.cpp +++ b/src/IO/mesh/impl/fortran/dim3/read_parameters.cpp @@ -8,223 +8,185 @@ specfem::mesh::parameters specfem::IO::mesh::impl::fortran::dim3::read_mesh_parameters( std::ifstream &stream, const specfem::MPI::MPI *mpi) { - // Initialize test parameter - int itest = -9999; - - // Read testparamater - specfem::IO::fortran_read_line(stream, &itest); + // Creating aliases for Array Reading functions + using specfem::IO::mesh::impl::fortran::dim3::check_read_test_value; + using specfem::IO::mesh::impl::fortran::dim3::try_read_line; - // Throw error if test parameter is not correctly read - if (itest != 9999) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - } + // Initialize test parameter + check_read_test_value(stream, 9999); - // Initialize flags - bool acoustic_simulation; - bool elastic_simulation; - bool poroelastic_simulation; - bool anisotropy; - bool stacey_abc; - bool pml_abc; - bool approximate_ocean_load; - bool use_mesh_coloring; + // Initialize parameters object + specfem::mesh::parameters parameters; // Read parameters - specfem::IO::fortran_read_line(stream, &acoustic_simulation); - specfem::IO::fortran_read_line(stream, &elastic_simulation); - specfem::IO::fortran_read_line(stream, &poroelastic_simulation); - specfem::IO::fortran_read_line(stream, &anisotropy); - specfem::IO::fortran_read_line(stream, &stacey_abc); - specfem::IO::fortran_read_line(stream, &pml_abc); - specfem::IO::fortran_read_line(stream, &approximate_ocean_load); - specfem::IO::fortran_read_line(stream, &use_mesh_coloring); + try_read_line("parameters_acoustic_simulation", stream, + ¶meters.acoustic_simulation); + try_read_line("parameters_elastic_simulation", stream, + ¶meters.elastic_simulation); + try_read_line("parameters_poroelastic_simulation", stream, + ¶meters.poroelastic_simulation); + try_read_line("parameters_anisotropy", stream, ¶meters.anisotropy); + try_read_line("parameters_stacey_abc", stream, ¶meters.stacey_abc); + try_read_line("parameters_pml_abc", stream, ¶meters.pml_abc); + try_read_line("parameters_approximate_ocean_load", stream, + ¶meters.approximate_ocean_load); + try_read_line("parameters_use_mesh_coloring", stream, + ¶meters.use_mesh_coloring); #ifndef NDEBUG // Print the read parameters - std::cout << "acoustic_simulation:...." << acoustic_simulation << std::endl; - std::cout << "elastic_simulation:....." << elastic_simulation << std::endl; - std::cout << "poroelastic_simulation:." << poroelastic_simulation + std::cout << "acoustic_simulation:...." << parameters.acoustic_simulation + << std::endl; + std::cout << "elastic_simulation:....." << parameters.elastic_simulation + << std::endl; + std::cout << "poroelastic_simulation:." << parameters.poroelastic_simulation + << std::endl; + std::cout << "anisotropy:............." << parameters.anisotropy << std::endl; + std::cout << "stacey_abc:............." << parameters.stacey_abc << std::endl; + std::cout << "pml_abc:................" << parameters.pml_abc << std::endl; + std::cout << "approximate_ocean_load:." << parameters.approximate_ocean_load << std::endl; - std::cout << "anisotropy:............." << anisotropy << std::endl; - std::cout << "stacey_abc:............." << stacey_abc << std::endl; - std::cout << "pml_abc:................" << pml_abc << std::endl; - std::cout << "approximate_ocean_load:." << approximate_ocean_load + std::cout << "use_mesh_coloring:......" << parameters.use_mesh_coloring << std::endl; - std::cout << "use_mesh_coloring:......" << use_mesh_coloring << std::endl; #endif // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); - - // Throw error if test parameter is not correctly read - if (itest != 9998) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; - - // Initialize first bout of integer parameters - int ndim; - int ngllx; - int nglly; - int ngllz; - int ngllsquare; - int nproc; + check_read_test_value(stream, 9998); // Read parameters - specfem::IO::fortran_read_line(stream, &ndim); - specfem::IO::fortran_read_line(stream, &ngllx); - specfem::IO::fortran_read_line(stream, &nglly); - specfem::IO::fortran_read_line(stream, &ngllz); - specfem::IO::fortran_read_line(stream, &ngllsquare); - specfem::IO::fortran_read_line(stream, &nproc); + try_read_line("mesh.parameters.ndim", stream, ¶meters.ndim); + try_read_line("mesh.parameters.ngllx", stream, ¶meters.ngllx); + try_read_line("mesh.parameters.nglly", stream, ¶meters.nglly); + try_read_line("mesh.parameters.ngllz", stream, ¶meters.ngllz); + try_read_line("mesh.parameters.ngllsquare", stream, ¶meters.ngllsquare); + try_read_line("mesh.parameters.nproc", stream, ¶meters.nproc); #ifndef NDEBUG // Print the read parameters - std::cout << "ndim:................." << ndim << std::endl; - std::cout << "ngllx:................" << ngllx << std::endl; - std::cout << "ngly:................." << nglly << std::endl; - std::cout << "ngllz:................" << ngllz << std::endl; - std::cout << "ngllsquare:..........." << ngllsquare << std::endl; - std::cout << "nproc:................" << nproc << std::endl; + std::cout << "ndim:................." << parameters.ndim << std::endl; + std::cout << "ngllx:................" << parameters.ngllx << std::endl; + std::cout << "ngly:................." << parameters.nglly << std::endl; + std::cout << "ngllz:................" << parameters.ngllz << std::endl; + std::cout << "ngllsquare:..........." << parameters.ngllsquare << std::endl; + std::cout << "nproc:................" << parameters.nproc << std::endl; #endif // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); - - // Throw error if test parameter is not correctly read - if (itest != 9997) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; - - // Initialize second bout of integer parameters - int nspec; - int nspec_poro; - int nglob; - int nglob_ocean; + check_read_test_value(stream, 9997); // Read parameters - specfem::IO::fortran_read_line(stream, &nspec); - specfem::IO::fortran_read_line(stream, &nspec_poro); - specfem::IO::fortran_read_line(stream, &nglob); - specfem::IO::fortran_read_line(stream, &nglob_ocean); + try_read_line("mesh.parameters.nspec", stream, ¶meters.nspec); + try_read_line("mesh.parameters.nspec_poro", stream, ¶meters.nspec_poro); + try_read_line("mesh.parameters.nglob", stream, ¶meters.nglob); + try_read_line("mesh.parameters.nglob_ocean", stream, ¶meters.nglob_ocean); #ifndef NDEBUG // Print the read parameters - std::cout << "nspec:................." << nspec << std::endl; - std::cout << "nspec_poro:............" << nspec_poro << std::endl; - std::cout << "nglob:................." << nglob << std::endl; - std::cout << "nglob_ocean:............" << nglob_ocean << std::endl; + std::cout << "nspec:................." << parameters.nspec << std::endl; + std::cout << "nspec_poro:............" << parameters.nspec_poro << std::endl; + std::cout << "nglob:................." << parameters.nglob << std::endl; + std::cout << "nglob_ocean:............" << parameters.nglob_ocean + << std::endl; #endif // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); + check_read_test_value(stream, 9996); - // Throw error if test parameter is not correctly read - if (itest != 9996) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; - - int nspec2D_bottom; - int nspec2D_top; - int nspec2D_xmin; - int nspec2D_xmax; - int nspec2D_ymin; - int nspec2D_ymax; - int nspec_irregular; - - specfem::IO::fortran_read_line(stream, &nspec2D_bottom); - specfem::IO::fortran_read_line(stream, &nspec2D_top); - specfem::IO::fortran_read_line(stream, &nspec2D_xmin); - specfem::IO::fortran_read_line(stream, &nspec2D_xmax); - specfem::IO::fortran_read_line(stream, &nspec2D_ymin); - specfem::IO::fortran_read_line(stream, &nspec2D_ymax); - specfem::IO::fortran_read_line(stream, &nspec_irregular); + // Read parameters + try_read_line("mesh.parameters.nspec2D_bottom", stream, + ¶meters.nspec2D_bottom); + try_read_line("mesh.parameters.nspec2D_top", stream, ¶meters.nspec2D_top); + try_read_line("mesh.parameters.nspec2D_xmin", stream, + ¶meters.nspec2D_xmin); + try_read_line("mesh.parameters.nspec2D_xmax", stream, + ¶meters.nspec2D_xmax); + try_read_line("mesh.parameters.nspec2D_ymin", stream, + ¶meters.nspec2D_ymin); + try_read_line("mesh.parameters.nspec2D_ymax", stream, + ¶meters.nspec2D_ymax); + try_read_line("mesh.parameters.nspec_irregular", stream, + ¶meters.nspec_irregular); #ifndef NDEBUG - std::cout << "nspec2D_bottom:........" << nspec2D_bottom << std::endl; - std::cout << "nspec2D_top:..........." << nspec2D_top << std::endl; - std::cout << "nspec2D_xmin:.........." << nspec2D_xmin << std::endl; - std::cout << "nspec2D_xmax:.........." << nspec2D_xmax << std::endl; - std::cout << "nspec2D_ymin:.........." << nspec2D_ymin << std::endl; - std::cout << "nspec2D_ymax:.........." << nspec2D_ymax << std::endl; - std::cout << "nspec_irregular:......." << nspec_irregular << std::endl; + std::cout << "nspec2D_bottom:........" << parameters.nspec2D_bottom + << std::endl; + std::cout << "nspec2D_top:..........." << parameters.nspec2D_top << std::endl; + std::cout << "nspec2D_xmin:.........." << parameters.nspec2D_xmin + << std::endl; + std::cout << "nspec2D_xmax:.........." << parameters.nspec2D_xmax + << std::endl; + std::cout << "nspec2D_ymin:.........." << parameters.nspec2D_ymin + << std::endl; + std::cout << "nspec2D_ymax:.........." << parameters.nspec2D_ymax + << std::endl; + std::cout << "nspec_irregular:......." << parameters.nspec_irregular + << std::endl; #endif // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); - - // Throw error if test parameter is not correctly read - if (itest != 9995) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; - - int num_neighbors; - int nfaces_surface; - int num_abs_boundary_faces; - int num_free_surface_faces; - int num_coupling_ac_el_faces; - int num_coupling_ac_po_faces; - int num_coupling_el_po_faces; - int num_coupling_po_el_faces; - int num_interfaces_ext_mesh; - int max_nibool_interfaces_ext_mesh; + check_read_test_value(stream, 9995); - specfem::IO::fortran_read_line(stream, &num_neighbors); - specfem::IO::fortran_read_line(stream, &nfaces_surface); - specfem::IO::fortran_read_line(stream, &num_abs_boundary_faces); - specfem::IO::fortran_read_line(stream, &num_free_surface_faces); - specfem::IO::fortran_read_line(stream, &num_coupling_ac_el_faces); - specfem::IO::fortran_read_line(stream, &num_coupling_ac_po_faces); - specfem::IO::fortran_read_line(stream, &num_coupling_el_po_faces); - specfem::IO::fortran_read_line(stream, &num_coupling_po_el_faces); - specfem::IO::fortran_read_line(stream, &num_interfaces_ext_mesh); - specfem::IO::fortran_read_line(stream, &max_nibool_interfaces_ext_mesh); + // Read parameters + try_read_line("mesh.parameters.num_neighbors", stream, + ¶meters.num_neighbors); + try_read_line("mesh.parameters.nfaces_surface", stream, + ¶meters.nfaces_surface); + try_read_line("mesh.parameters.num_abs_boundary_faces", stream, + ¶meters.num_abs_boundary_faces); + try_read_line("mesh.parameters.num_free_surface_faces", stream, + ¶meters.num_free_surface_faces); + try_read_line("mesh.parameters.num_coupling_ac_el_faces", stream, + ¶meters.num_coupling_ac_el_faces); + try_read_line("mesh.parameters.num_coupling_ac_po_faces", stream, + ¶meters.num_coupling_ac_po_faces); + try_read_line("mesh.parameters.num_coupling_el_po_faces", stream, + ¶meters.num_coupling_el_po_faces); + try_read_line("mesh.parameters.num_coupling_po_el_faces", stream, + ¶meters.num_coupling_po_el_faces); + try_read_line("mesh.parameters.num_interfaces_ext_mesh", stream, + ¶meters.num_interfaces_ext_mesh); + try_read_line("mesh.parameters.max_nibool_interfaces_ext_mesh", stream, + ¶meters.max_nibool_interfaces_ext_mesh); #ifndef NDEBUG - std::cout << "num_neighbors:.........." << num_neighbors << std::endl; - std::cout << "nfaces_surface:........." << nfaces_surface << std::endl; - std::cout << "num_abs_boundary_faces:." << num_abs_boundary_faces - << std::endl; - std::cout << "num_free_surface_faces:." << num_free_surface_faces - << std::endl; - std::cout << "num_coupling_ac_el_faces:" << num_coupling_ac_el_faces + std::cout << "num_neighbors:.........." << parameters.num_neighbors << std::endl; - std::cout << "num_coupling_ac_po_faces:" << num_coupling_ac_po_faces + std::cout << "nfaces_surface:........." << parameters.nfaces_surface << std::endl; - std::cout << "num_coupling_el_po_faces:" << num_coupling_el_po_faces + std::cout << "num_abs_boundary_faces:." << parameters.num_abs_boundary_faces << std::endl; - std::cout << "num_coupling_po_el_faces:" << num_coupling_po_el_faces + std::cout << "num_free_surface_faces:." << parameters.num_free_surface_faces << std::endl; - std::cout << "num_interfaces_ext_mesh:." << num_interfaces_ext_mesh + std::cout << "num_coupling_ac_el_faces:" + << parameters.num_coupling_ac_el_faces << std::endl; + std::cout << "num_coupling_ac_po_faces:" + << parameters.num_coupling_ac_po_faces << std::endl; + std::cout << "num_coupling_el_po_faces:" + << parameters.num_coupling_el_po_faces << std::endl; + std::cout << "num_coupling_po_el_faces:" + << parameters.num_coupling_po_el_faces << std::endl; + std::cout << "num_interfaces_ext_mesh:." << parameters.num_interfaces_ext_mesh << std::endl; std::cout << "max_nibool_interfaces_ext_mesh:" - << max_nibool_interfaces_ext_mesh << std::endl; + << parameters.max_nibool_interfaces_ext_mesh << std::endl; #endif // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); + check_read_test_value(stream, 9994); - // Throw error if test parameter is not correctly read - if (itest != 9994) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; + // Read parameters + try_read_line("mesh.parameters.nspec_inner_acoustic", stream, + ¶meters.nspec_inner_acoustic); + try_read_line("mesh.parameters.nspec_outer_acoustic", stream, + ¶meters.nspec_outer_acoustic); + try_read_line("mesh.parameters.nspec_inner_elastic", stream, + ¶meters.nspec_inner_elastic); + try_read_line("mesh.parameters.nspec_outer_elastic", stream, + ¶meters.nspec_outer_elastic); + try_read_line("mesh.parameters.nspec_inner_poroelastic", stream, + ¶meters.nspec_inner_poroelastic); + try_read_line("mesh.parameters.nspec_outer_poroelastic", stream, + ¶meters.nspec_outer_poroelastic); int nspec_inner_acoustic; int nspec_outer_acoustic; @@ -233,127 +195,61 @@ specfem::IO::mesh::impl::fortran::dim3::read_mesh_parameters( int nspec_inner_poroelastic; int nspec_outer_poroelastic; - specfem::IO::fortran_read_line(stream, &nspec_inner_acoustic); - specfem::IO::fortran_read_line(stream, &nspec_outer_acoustic); - specfem::IO::fortran_read_line(stream, &nspec_inner_elastic); - specfem::IO::fortran_read_line(stream, &nspec_outer_elastic); - specfem::IO::fortran_read_line(stream, &nspec_inner_poroelastic); - specfem::IO::fortran_read_line(stream, &nspec_outer_poroelastic); - #ifndef NDEBUG - std::cout << "nspec_inner_acoustic:..." << nspec_inner_acoustic << std::endl; - std::cout << "nspec_outer_acoustic:..." << nspec_outer_acoustic << std::endl; - std::cout << "nspec_inner_elastic:...." << nspec_inner_elastic << std::endl; - std::cout << "nspec_outer_elastic:...." << nspec_outer_elastic << std::endl; - std::cout << "nspec_inner_poroelastic:." << nspec_inner_poroelastic + std::cout << "nspec_inner_acoustic:..." << parameters.nspec_inner_acoustic << std::endl; - std::cout << "nspec_outer_poroelastic:." << nspec_outer_poroelastic + std::cout << "nspec_outer_acoustic:..." << parameters.nspec_outer_acoustic << std::endl; -#endif - - // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); - - // Throw error if test parameter is not correctly read - if (itest != 9993) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; - - int num_phase_ispec_acoustic; - int num_phase_ispec_elastic; - int num_phase_ispec_poroelastic; - int num_colors_inner_acoustic; - int num_colors_outer_acoustic; - int num_colors_inner_elastic; - int num_colors_outer_elastic; - - specfem::IO::fortran_read_line(stream, &num_phase_ispec_acoustic); - specfem::IO::fortran_read_line(stream, &num_phase_ispec_elastic); - specfem::IO::fortran_read_line(stream, &num_phase_ispec_poroelastic); - specfem::IO::fortran_read_line(stream, &num_colors_inner_acoustic); - specfem::IO::fortran_read_line(stream, &num_colors_outer_acoustic); - specfem::IO::fortran_read_line(stream, &num_colors_inner_elastic); - specfem::IO::fortran_read_line(stream, &num_colors_outer_elastic); - -#ifndef NDEBUG - std::cout << "num_phase_ispec_acoustic:" << num_phase_ispec_acoustic - << std::endl; - std::cout << "num_phase_ispec_elastic:." << num_phase_ispec_elastic - << std::endl; - std::cout << "num_phase_ispec_poroelastic:" << num_phase_ispec_poroelastic - << std::endl; - std::cout << "num_colors_inner_acoustic:" << num_colors_inner_acoustic + std::cout << "nspec_inner_elastic:...." << parameters.nspec_inner_elastic << std::endl; - std::cout << "num_colors_outer_acoustic:" << num_colors_outer_acoustic + std::cout << "nspec_outer_elastic:...." << parameters.nspec_outer_elastic << std::endl; - std::cout << "num_colors_inner_elastic:" << num_colors_inner_elastic + std::cout << "nspec_inner_poroelastic:." << parameters.nspec_inner_poroelastic << std::endl; - std::cout << "num_colors_outer_elastic:" << num_colors_outer_elastic + std::cout << "nspec_outer_poroelastic:." << parameters.nspec_outer_poroelastic << std::endl; #endif // Read test parameter - specfem::IO::fortran_read_line(stream, &itest); + check_read_test_value(stream, 9993); + + // Read parameters + try_read_line("mesh.parameters.num_phase_ispec_acoustic", stream, + ¶meters.num_phase_ispec_acoustic); + try_read_line("mesh.parameters.num_phase_ispec_elastic", stream, + ¶meters.num_phase_ispec_elastic); + try_read_line("mesh.parameters.num_phase_ispec_poroelastic", stream, + ¶meters.num_phase_ispec_poroelastic); + try_read_line("mesh.parameters.num_colors_inner_acoustic", stream, + ¶meters.num_colors_inner_acoustic); + try_read_line("mesh.parameters.num_colors_outer_acoustic", stream, + ¶meters.num_colors_outer_acoustic); + try_read_line("mesh.parameters.num_colors_inner_elastic", stream, + ¶meters.num_colors_inner_elastic); + try_read_line("mesh.parameters.num_colors_outer_elastic", stream, + ¶meters.num_colors_outer_elastic); + +#ifndef NDEBUG + std::cout << "num_phase_ispec_acoustic:" + << parameters.num_phase_ispec_acoustic << std::endl; + std::cout << "num_phase_ispec_elastic:." << parameters.num_phase_ispec_elastic + << std::endl; + std::cout << "num_phase_ispec_poroelastic:" + << parameters.num_phase_ispec_poroelastic << std::endl; + std::cout << "num_colors_inner_acoustic:" + << parameters.num_colors_inner_acoustic << std::endl; + std::cout << "num_colors_outer_acoustic:" + << parameters.num_colors_outer_acoustic << std::endl; + std::cout << "num_colors_inner_elastic:" + << parameters.num_colors_inner_elastic << std::endl; + std::cout << "num_colors_outer_elastic:" + << parameters.num_colors_outer_elastic << std::endl; +#endif - // Throw error if test parameter is not correctly read - if (itest != 9992) { - std::ostringstream error_message; - error_message << "Error reading mesh parameters: " << itest - << " // FILE:LINE " << __FILE__ << ":" << __LINE__ << "\n"; - throw std::runtime_error(error_message.str()); - }; + /// Read test parameter + check_read_test_value(stream, 9992); mpi->sync_all(); - return { acoustic_simulation, - elastic_simulation, - poroelastic_simulation, - anisotropy, - stacey_abc, - pml_abc, - approximate_ocean_load, - use_mesh_coloring, - ndim, - ngllx, - nglly, - ngllz, - ngllsquare, - nproc, - nspec, - nspec_poro, - nglob, - nglob_ocean, - nspec2D_bottom, - nspec2D_top, - nspec2D_xmin, - nspec2D_xmax, - nspec2D_ymin, - nspec2D_ymax, - nspec_irregular, - num_neighbors, - nfaces_surface, - num_abs_boundary_faces, - num_free_surface_faces, - num_coupling_ac_el_faces, - num_coupling_ac_po_faces, - num_coupling_el_po_faces, - num_coupling_po_el_faces, - num_interfaces_ext_mesh, - max_nibool_interfaces_ext_mesh, - nspec_inner_acoustic, - nspec_outer_acoustic, - nspec_inner_elastic, - nspec_outer_elastic, - nspec_inner_poroelastic, - nspec_outer_poroelastic, - num_phase_ispec_acoustic, - num_phase_ispec_elastic, - num_phase_ispec_poroelastic, - num_colors_inner_acoustic, - num_colors_outer_acoustic, - num_colors_inner_elastic, - num_colors_outer_elastic }; + return parameters; } diff --git a/src/IO/mesh/impl/fortran/dim3/read_partial_derivatives.cpp b/src/IO/mesh/impl/fortran/dim3/read_partial_derivatives.cpp index d1bfa0eb3..8f75a96bd 100644 --- a/src/IO/mesh/impl/fortran/dim3/read_partial_derivatives.cpp +++ b/src/IO/mesh/impl/fortran/dim3/read_partial_derivatives.cpp @@ -35,8 +35,8 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( // Read all elements at once try { - specfem::IO::mesh::impl::fortran::dim3::read_array( - stream, partial_derivatives.xix); + specfem::IO::mesh::impl::fortran::dim3::read_array(stream, + partial_derivatives.xix); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -46,8 +46,8 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( - stream, partial_derivatives.xiy); + specfem::IO::mesh::impl::fortran::dim3::read_array(stream, + partial_derivatives.xiy); } catch (std::runtime_error &e) { std::ostringstream error_message; error_message << "Error reading xiy from database file:\n" @@ -56,8 +56,8 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( - stream, partial_derivatives.xiz); + specfem::IO::mesh::impl::fortran::dim3::read_array(stream, + partial_derivatives.xiz); } catch (std::runtime_error &e) { std::ostringstream error_message; error_message << "Error reading xiz from database file:\n" @@ -66,7 +66,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( + specfem::IO::mesh::impl::fortran::dim3::read_array( stream, partial_derivatives.etax); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -76,7 +76,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( + specfem::IO::mesh::impl::fortran::dim3::read_array( stream, partial_derivatives.etay); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -86,7 +86,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( + specfem::IO::mesh::impl::fortran::dim3::read_array( stream, partial_derivatives.etaz); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -96,7 +96,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( + specfem::IO::mesh::impl::fortran::dim3::read_array( stream, partial_derivatives.gammax); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -106,7 +106,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( + specfem::IO::mesh::impl::fortran::dim3::read_array( stream, partial_derivatives.gammay); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -116,7 +116,7 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( } try { - specfem::IO::mesh::impl::fortran::dim3::read_array( + specfem::IO::mesh::impl::fortran::dim3::read_array( stream, partial_derivatives.gammaz); } catch (std::runtime_error &e) { std::ostringstream error_message; @@ -124,4 +124,14 @@ void specfem::IO::mesh::impl::fortran::dim3::read_partial_derivatives( << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; throw std::runtime_error(error_message.str()); } + + try { + specfem::IO::mesh::impl::fortran::dim3::read_array( + stream, partial_derivatives.jacobian); + } catch (std::runtime_error &e) { + std::ostringstream error_message; + error_message << "Error reading jacobian from database file:\n" + << e.what() << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } } diff --git a/src/IO/mesh/impl/fortran/dim3/utilities.cpp b/src/IO/mesh/impl/fortran/dim3/utilities.cpp new file mode 100644 index 000000000..88139a2a3 --- /dev/null +++ b/src/IO/mesh/impl/fortran/dim3/utilities.cpp @@ -0,0 +1,32 @@ +#include "IO/mesh/impl/fortran/dim3/interface.hpp" +#include +#include +#include +#include + +// Rewrite check_read_test_value to use try_read_line +void specfem::IO::mesh::impl::fortran::dim3::check_read_test_value( + std::ifstream &stream, int test_value) { + // Read test value that should be value + int value; + try_read_line("check_read_test_value", stream, &value); + if (test_value != value) { + std::ostringstream error_message; + error_message << "Test value (" << test_value << ") != read value (" + << value << "). " + << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } +} + +void specfem::IO::mesh::impl::fortran::dim3::check_values(std::string message, + int value, + int expected) { + if (value != expected) { + std::ostringstream error_message; + error_message << message << " value (" << value << ") != expected value (" + << expected << "). " + << "(" << __FILE__ << ":" << __LINE__ << ")"; + throw std::runtime_error(error_message.str()); + } +} diff --git a/src/compute/compute_properties.cpp b/src/compute/compute_properties.cpp index f21b2c0f1..5af22bdfe 100644 --- a/src/compute/compute_properties.cpp +++ b/src/compute/compute_properties.cpp @@ -1,9 +1,11 @@ #include "compute/properties/properties.hpp" +#include "enumerations/dimension.hpp" specfem::compute::properties::properties( const int nspec, const int ngllz, const int ngllx, const specfem::compute::element_types &element_types, - const specfem::mesh::materials &materials, const bool has_gll_model) { + const specfem::mesh::materials &materials, + const bool has_gll_model) { this->nspec = nspec; this->ngllz = ngllz; diff --git a/src/mesh/dim2/coupled_interfaces/coupled_interfaces.cpp b/src/mesh/dim2/coupled_interfaces/coupled_interfaces.cpp index 4e965b904..9e10afb79 100644 --- a/src/mesh/dim2/coupled_interfaces/coupled_interfaces.cpp +++ b/src/mesh/dim2/coupled_interfaces/coupled_interfaces.cpp @@ -1,6 +1,7 @@ #include "mesh/dim2/coupled_interfaces/coupled_interfaces.hpp" #include "mesh/dim2/coupled_interfaces/interface_container.hpp" #include "mesh/dim2/coupled_interfaces/interface_container.tpp" +#include "mesh/mesh_base.hpp" // specfem::mesh::coupled_interfaces::coupled_interfaces( // specfem::mesh::interface_container // elastic_poroelastic) {} -template template -std::variant, - specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::acoustic, - specfem::element::medium_tag::poroelastic>, - specfem::mesh::interface_container< - DimensionType, specfem::element::medium_tag::elastic, - specfem::element::medium_tag::poroelastic> > -specfem::mesh::coupled_interfaces::coupled_interfaces::get() - const { +std::variant< + specfem::mesh::interface_container, + specfem::mesh::interface_container< + specfem::dimension::type::dim2, specfem::element::medium_tag::acoustic, + specfem::element::medium_tag::poroelastic>, + specfem::mesh::interface_container< + specfem::dimension::type::dim2, specfem::element::medium_tag::elastic, + specfem::element::medium_tag::poroelastic> > +specfem::mesh::coupled_interfaces< + specfem::dimension::type::dim2>::coupled_interfaces::get() const { if constexpr (Medium1 == specfem::element::medium_tag::elastic && Medium2 == specfem::element::medium_tag::acoustic) { return elastic_acoustic; diff --git a/src/mesh/dim2/coupled_interfaces/interface_container.cpp b/src/mesh/dim2/coupled_interfaces/interface_container.cpp index eb83cff04..05a03f414 100644 --- a/src/mesh/dim2/coupled_interfaces/interface_container.cpp +++ b/src/mesh/dim2/coupled_interfaces/interface_container.cpp @@ -1,5 +1,6 @@ #include "mesh/dim2/coupled_interfaces/interface_container.hpp" #include "mesh/dim2/coupled_interfaces/interface_container.tpp" +#include "mesh/mesh_base.hpp" #define GENERATE_INTERFACE_CONTAINER(DimensionType) \ template class specfem::mesh::interface_container< \ DimensionType, specfem::element::medium_tag::elastic, \ diff --git a/src/mesh/dim2/materials/materials.cpp b/src/mesh/dim2/materials/materials.cpp index 15ce0b3b5..dfda29fce 100644 --- a/src/mesh/dim2/materials/materials.cpp +++ b/src/mesh/dim2/materials/materials.cpp @@ -11,7 +11,8 @@ std::variant< specfem::element::property_tag::anisotropic>, specfem::medium::material > -specfem::mesh::materials::operator[](const int index) const { +specfem::mesh::materials::operator[]( + const int index) const { const auto &material_specification = this->material_index_mapping(index); diff --git a/src/mesh/dim2/tags/tags.cpp b/src/mesh/dim2/tags/tags.cpp index 6f68d2967..7905b2873 100644 --- a/src/mesh/dim2/tags/tags.cpp +++ b/src/mesh/dim2/tags/tags.cpp @@ -2,7 +2,7 @@ #include "enumerations/dimension.hpp" specfem::mesh::tags::tags( - const specfem::mesh::materials &materials, + const specfem::mesh::materials &materials, const specfem::mesh::boundaries &boundaries) { diff --git a/src/mesh/dim3/boundaries/absorbing_boundary.cpp b/src/mesh/dim3/boundaries/absorbing_boundary.cpp new file mode 100644 index 000000000..d0ffe6039 --- /dev/null +++ b/src/mesh/dim3/boundaries/absorbing_boundary.cpp @@ -0,0 +1,41 @@ +#include "mesh/mesh.hpp" +#include + +std::string +specfem::mesh::absorbing_boundary::print() + const { + + std::ostringstream message; + + // Print variables and metadata + message << "Absorbing Boundary Metadata:" + << "\n"; + message << "===============================================" + << "\n"; + message << " nelements:.............. " << nelements << "\n"; + message << " ngllsquare:............. " << ngllsquare << "\n"; + message << " num_abs_boundary_faces:. " << num_abs_boundary_faces << "\n"; + message << " acoustic:............... " << acoustic << "\n"; + message << " elastic:................ " << elastic << "\n"; + message << " nspec2D_xmin:........... " << nspec2D_xmin << "\n"; + message << " nspec2D_xmax:........... " << nspec2D_xmax << "\n"; + message << " nspec2D_ymin:........... " << nspec2D_ymin << "\n"; + message << " nspec2D_ymax:........... " << nspec2D_ymax << "\n"; + message << " NSPEC2D_BOTTOM:......... " << NSPEC2D_BOTTOM << "\n"; + message << " NSPEC2D_TOP:............ " << NSPEC2D_TOP << "\n"; + + // Print the absorbing ispec metadata + message << " Array sizes:" + << "\n"; + message << " -----------------------------------------------" + << "\n"; + message << " ispec:.................. " << ispec.extent(0) << "\n"; + message << " ijk:.................... " << ijk.extent(0) << " " + << ijk.extent(1) << " " << ijk.extent(2) << "\n"; + message << " jacobian2Dw:............ " << jacobian2Dw.extent(0) << " " + << jacobian2Dw.extent(1) << "\n"; + message << " normal:................. " << normal.extent(0) << " " + << normal.extent(1) << " " << normal.extent(2) << "\n"; + + return message.str(); +} diff --git a/src/mesh/dim3/boundaries/free_surface.cpp b/src/mesh/dim3/boundaries/free_surface.cpp new file mode 100644 index 000000000..678b7bf88 --- /dev/null +++ b/src/mesh/dim3/boundaries/free_surface.cpp @@ -0,0 +1,33 @@ +#include "mesh/mesh.hpp" +#include +#include + +std::string +specfem::mesh::free_surface::print() const { + + std::ostringstream message; + + // Print variables and metadata + message << "Absorbing Boundary Metadata:" + << "\n"; + message << "================================================" + << "\n"; + message << " nelements:.............. " << nelements << "\n"; + message << " ngllsquare:............. " << ngllsquare << "\n"; + message << " num_free_surface_faces:. " << num_free_surface_faces << "\n"; + + // Print the absorbing ispec metadata + message << " Array sizes:" + << "\n"; + message << " -----------------------------------------------" + << "\n"; + message << " ispec:.................. " << ispec.extent(0) << "\n"; + message << " ijk:.................... " << ijk.extent(0) << " " + << ijk.extent(1) << " " << ijk.extent(2) << "\n"; + message << " jacobian2Dw:............ " << jacobian2Dw.extent(0) << " " + << jacobian2Dw.extent(1) << "\n"; + message << " normal:................. " << normal.extent(0) << " " + << normal.extent(1) << " " << normal.extent(2) << "\n"; + + return message.str(); +} diff --git a/src/mesh/dim3/coordinates/coordinates.cpp b/src/mesh/dim3/coordinates/coordinates.cpp index bae3be85e..93ac1efcb 100644 --- a/src/mesh/dim3/coordinates/coordinates.cpp +++ b/src/mesh/dim3/coordinates/coordinates.cpp @@ -1,23 +1,31 @@ #include "mesh/dim3/coordinates/coordinates.hpp" #include -void specfem::mesh::coordinates::print() const { - std::cout << "Coordinates parameters:\n" - << "------------------------------\n" - << "Number of spectral elements: " << nspec << "\n" - << "Number of global nodes:..... " << nglob << "\n" - << "Number of GLLX:............. " << ngllx << "\n" - << "Number of GLLY:............. " << nglly << "\n" - << "Number of GLLZ:............. " << ngllz << "\n" - << "------------------------------\n"; +std::string +specfem::mesh::coordinates::print() const { + + std::ostringstream message; + message << "Coordinates parameters:\n" + << "------------------------------\n" + << "Number of spectral elements: " << nspec << "\n" + << "Number of global nodes:..... " << nglob << "\n" + << "Number of GLLX:............. " << ngllx << "\n" + << "Number of GLLY:............. " << nglly << "\n" + << "Number of GLLZ:............. " << ngllz << "\n" + << "------------------------------\n"; + + return message.str(); } -void specfem::mesh::coordinates::print( +std::string specfem::mesh::coordinates::print( int iglob) const { - std::cout << "Coordinates parameters for global node " << iglob << ":\n" - << "------------------------------\n" - << "x: " << x(iglob) << "\n" - << "y: " << y(iglob) << "\n" - << "z: " << z(iglob) << "\n" - << "------------------------------\n"; + std::ostringstream message; + message << "Coordinates parameters for global node " << iglob << ":\n" + << "------------------------------\n" + << "x: " << x(iglob) << "\n" + << "y: " << y(iglob) << "\n" + << "z: " << z(iglob) << "\n" + << "------------------------------\n"; + + return message.str(); } diff --git a/src/mesh/dim3/coupled_interfaces/coupled_interfaces.cpp b/src/mesh/dim3/coupled_interfaces/coupled_interfaces.cpp new file mode 100644 index 000000000..d07a1ca15 --- /dev/null +++ b/src/mesh/dim3/coupled_interfaces/coupled_interfaces.cpp @@ -0,0 +1,110 @@ +#include "mesh/mesh.hpp" +#include + +std::string +specfem::mesh::coupled_interfaces::print() + const { + + std::ostringstream message; + + if (acoustic_elastic) { + message << "Acoustic Elastic Interface Metadata:" + << "\n"; + message << "================================================" + << "\n"; + message << " ngllsquare:............... " << ngllsquare << "\n"; + message << " num_coupling_ac_el_faces:. " << num_coupling_ac_el_faces + << "\n"; + + // Print the acoustic elastic interface metadata + message << " Array sizes:" + << "\n"; + message << " -----------------------------------------------" + << "\n"; + message << " ispec:.................. " + << acoustic_elastic_interface.ispec.extent(0) << "\n"; + message << " ijk:.................... " + << acoustic_elastic_interface.ijk.extent(0) << " " + << acoustic_elastic_interface.ijk.extent(1) << " " + << acoustic_elastic_interface.ijk.extent(2) << "\n"; + message << " jacobian2Dw:............ " + << acoustic_elastic_interface.jacobian2Dw.extent(0) << " " + << acoustic_elastic_interface.jacobian2Dw.extent(1) << "\n"; + message << " normal:................. " + << acoustic_elastic_interface.normal.extent(0) << " " + << acoustic_elastic_interface.normal.extent(1) << " " + << acoustic_elastic_interface.normal.extent(2) << "\n"; + } else { + message << "No acoustic elastic interfaces" + << "\n"; + } + + if (acoustic_poroelastic) { + + message << "Acoustic Poroelastic Interface Metadata:" + << "\n"; + message << "================================================" + << "\n"; + message << " ngllsquare:............... " << ngllsquare << "\n"; + message << " num_coupling_ac_po_faces:. " << num_coupling_ac_po_faces + << "\n"; + + // Print the acoustic poroelastic interface metadata + message << " Array sizes:" + << "\n"; + message << " -----------------------------------------------" + << "\n"; + message << " ispec:.................. " + << acoustic_poroelastic_interface.ispec.extent(0) << "\n"; + message << " ijk:.................... " + << acoustic_poroelastic_interface.ijk.extent(0) << " " + << acoustic_poroelastic_interface.ijk.extent(1) << " " + << acoustic_poroelastic_interface.ijk.extent(2) << "\n"; + message << " jacobian2Dw:............ " + << acoustic_poroelastic_interface.jacobian2Dw.extent(0) << " " + << acoustic_poroelastic_interface.jacobian2Dw.extent(1) << "\n"; + message << " normal:................. " + << acoustic_poroelastic_interface.normal.extent(0) << " " + << acoustic_poroelastic_interface.normal.extent(1) << " " + << acoustic_poroelastic_interface.normal.extent(2) << "\n"; + } else { + message << "No acoustic poroelastic interfaces" + << "\n"; + } + + if (elastic_poroelastic) { + + message << "Elastic Poroelastic Interface Metadata:" + << "\n"; + message << "================================================" + << "\n"; + message << " ngllsquare:............... " << ngllsquare << "\n"; + message << " num_coupling_el_po_faces:. " << num_coupling_el_po_faces + << "\n"; + + // Print the elastic poroelastic interface metadata + message << " Array sizes:" + << "\n"; + message << " -----------------------------------------------" + << "\n"; + message << " ispec:.................. " + << elastic_poroelastic_interface.ispec.extent(0) << "\n"; + message << " ijk:.................... " + << elastic_poroelastic_interface.ijk.extent(0) << " " + << elastic_poroelastic_interface.ijk.extent(1) << " " + << elastic_poroelastic_interface.ijk.extent(2) << "\n"; + message << " jacobian2Dw:............ " + << elastic_poroelastic_interface.jacobian2Dw.extent(0) << " " + << elastic_poroelastic_interface.jacobian2Dw.extent(1) << "\n"; + message << " normal:................. " + << elastic_poroelastic_interface.normal.extent(0) << " " + << elastic_poroelastic_interface.normal.extent(1) << " " + << elastic_poroelastic_interface.normal.extent(2) << "\n"; + + } else { + message << "No elastic poroelastic interfaces" + << "\n"; + } + + return message.str(); +} diff --git a/src/mesh/dim3/element_types/element_types.cpp b/src/mesh/dim3/element_types/element_types.cpp new file mode 100644 index 000000000..f75c70cbe --- /dev/null +++ b/src/mesh/dim3/element_types/element_types.cpp @@ -0,0 +1,134 @@ +#include "mesh/dim3/element_types/element_types.hpp" +#include "enumerations/dimension.hpp" +#include "enumerations/medium.hpp" +#include + +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; + + // Initialize the number of elements + this->nelastic = 0; + this->nacoustic = 0; + this->nporoelastic = 0; + + // Initialize the vectors + for (int ispec = 0; ispec < nspec; 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_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; + message << "Error: Spectral element " << ispec + << " is not assigned to any material"; + throw std::runtime_error(message.str()); + } + }; +} + +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(); +} + +std::string specfem::mesh::element_types::print( + const int ispec) const { + + 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 +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; + message << "Error: Elastic element " << i + << " does not exist.\nNumber " + "of elastic elements is :" + << nelastic << ". " + << "(" << __FILE__ << ":" << __LINE__ << ")\n"; + throw std::runtime_error(message.str()); + } + message << "Elastic element " << i << " is global element " + << ispec_elastic(i) << ".\n"; + } else if (MediumTag == specfem::element::medium_tag::acoustic) { + + if (i >= nacoustic) { + std::ostringstream message; + message << "Error: Acoustic element " << i + << " does not exist.\nNumber " + " of acoustic elements: " + << nacoustic << ". " + << "(" << __FILE__ << ":" << __LINE__ << ")\n"; + throw std::runtime_error(message.str()); + } + + message << "Acoustic element " << i << " is global element " + << ispec_acoustic(i) << ".\n"; + + } else if (MediumTag == specfem::element::medium_tag::poroelastic) { + + if (i >= nporoelastic) { + std::ostringstream message; + message << "Error: Poroelastic element " << i + << " does not exist.\nNumber " + " of poroelastic elements: " + << nporoelastic << ". " + << "(" << __FILE__ << ":" << __LINE__ << ")\n"; + throw std::runtime_error(message.str()); + } + + message << "Poroelastic element " << i << " is global element " + << ispec_poroelastic(i) << ".\n"; + } + + return message.str(); +} + +// Explicit instantiations +template std::string +specfem::mesh::element_types::print< + specfem::element::medium_tag::elastic>(const int) const; +template std::string +specfem::mesh::element_types::print< + specfem::element::medium_tag::acoustic>(const int) const; +template std::string +specfem::mesh::element_types::print< + specfem::element::medium_tag::poroelastic>(const int) const; diff --git a/src/mesh/dim3/mapping/mapping.cpp b/src/mesh/dim3/mapping/mapping.cpp index af020f3e6..7661cc341 100644 --- a/src/mesh/dim3/mapping/mapping.cpp +++ b/src/mesh/dim3/mapping/mapping.cpp @@ -1,33 +1,55 @@ #include "mesh/dim3/mapping/mapping.hpp" #include +#include -void specfem::mesh::mapping::print() const { - std::cout << "Mapping parameters:\n" - << "------------------------------\n" - << "Number of spectral elements: " << nspec << "\n" - << "Number of global nodes: " << nglob << "\n" - << "Number of irregular spectral elements: " << nspec_irregular - << "\n" - << "Number of GLLX: " << ngllx << "\n" - << "Number of GLLY: " << nglly << "\n" - << "Number of GLLZ: " << ngllz << "\n" - << "------------------------------\n"; +std::string +specfem::mesh::mapping::print() const { + + std::ostringstream message; + + message << "Mapping parameters:\n" + << "------------------------------\n" + << "Number of spectral elements: " << nspec << "\n" + << "Number of global nodes: " << nglob << "\n" + << "Number of irregular spectral elements: " << nspec_irregular + << "\n" + << "Number of GLLX: " << ngllx << "\n" + << "Number of GLLY: " << nglly << "\n" + << "Number of GLLZ: " << ngllz << "\n" + << "------------------------------\n"; + + return message.str(); } -void specfem::mesh::mapping::print( - int ispec) const { - std::cout << "Mapping parameters for spectral element " << ispec << ":\n" - << "------------------------------\n" - << "ibool:\n"; +std::string +specfem::mesh::mapping::print(int ispec) const { + + std::ostringstream message; + + message << "Mapping parameters for spectral element " << ispec << ":\n" + << "--------------------------------------------------\n" + << "\n" + << " |---> igllx\n" + << " |\n" + << " V\n" + << "iglly\n" + << "\n" + << "ibool:\n"; for (int igllz = 0; igllz < ngllz; igllz++) { + message << "igllz=" << igllz << ": "; for (int iglly = 0; iglly < nglly; iglly++) { + if (iglly > 0) { + message << " "; + } for (int igllx = 0; igllx < ngllx; igllx++) { - std::cout << ibool(ispec, igllx, iglly, igllz) << " "; + message << ibool(ispec, igllx, iglly, igllz) << " "; } - std::cout << "\n"; + message << "\n"; } - std::cout << "\n"; + message << "\n"; } - std::cout << "------------------------------\n"; + message << "------------------------------\n"; + + return message.str(); } diff --git a/src/mesh/dim3/materials/materials.cpp b/src/mesh/dim3/materials/materials.cpp new file mode 100644 index 000000000..3368ca9df --- /dev/null +++ b/src/mesh/dim3/materials/materials.cpp @@ -0,0 +1,80 @@ +#include "mesh/dim3/materials/materials.hpp" +#include "enumerations/dimension.hpp" +#include "kokkos_abstractions.h" +#include "specfem_setup.hpp" +#include +#include +#include +#include + +std::string specfem::mesh::materials::print() { + + std::ostringstream message; + message << "---------------------------------" + << "\n"; + message << "Materials: " + << "\n"; + message << "---------------------------------" + << "\n"; + message << "nspec: " << nspec << "\n"; + message << "ngllx: " << ngllx << "\n"; + message << "nglly: " << nglly << "\n"; + message << "ngllz: " << ngllz << "\n"; + message << "acoustic: " << acoustic << "\n"; + message << "elastic: " << elastic << "\n"; + message << "poroelastic: " << poroelastic << "\n"; + message << "anisotropic: " << anisotropic << "\n"; + + // Initialize min and max with appropriate extreme values + type_real kappa_min = std::numeric_limits::max(); + type_real kappa_max = std::numeric_limits::lowest(); + + Kokkos::parallel_reduce( + "kappa_min_max", + specfem::kokkos::HostRange(0, nspec * ngllx * nglly * ngllz), + KOKKOS_LAMBDA(const int &i, type_real &lmin, type_real &lmax) { + // Compute 3D indices from 1D index i + const int ispec = i / (ngllx * nglly * ngllz); + const int j = (i / (ngllx * ngllz)) % nglly; + const int k = (i / ngllx) % ngllz; + const int l = i % ngllx; + + const type_real val = kappa(ispec, j, k, l); + + // Update local min and max + lmin = val < lmin ? val : lmin; + lmax = val > lmax ? val : lmax; + }, + Kokkos::Min(kappa_min), Kokkos::Max(kappa_max)); + + message << "Kappa min/max: " << kappa_min << "/" << kappa_max << "\n"; + + // Initialize min and max with appropriate extreme values + type_real mu_min = std::numeric_limits::max(); + type_real mu_max = std::numeric_limits::lowest(); + + Kokkos::parallel_reduce( + "mu_min_max", + specfem::kokkos::HostRange(0, nspec * ngllx * nglly * ngllz), + KOKKOS_LAMBDA(const int &i, type_real &lmin, type_real &lmax) { + // Compute 3D indices from 1D index i + const int ispec = i / (ngllx * nglly * ngllz); + const int j = (i / (ngllx * ngllz)) % nglly; + const int k = (i / ngllx) % ngllz; + const int l = i % ngllx; + + const type_real val = mu(ispec, j, k, l); + + // Update local min and max + lmin = val < lmin ? val : lmin; + lmax = val > lmax ? val : lmax; + }, + Kokkos::Min(mu_min), Kokkos::Max(mu_max)); + + message << "Mu min/max: " << mu_min << "/" << mu_max << "\n"; + + message << "---------------------------------" + << "\n"; + + return message.str(); +} diff --git a/src/mesh/dim3/mesh.cpp b/src/mesh/dim3/mesh.cpp index c4c88c395..b75380089 100644 --- a/src/mesh/dim3/mesh.cpp +++ b/src/mesh/dim3/mesh.cpp @@ -10,35 +10,37 @@ #include #include -std::string specfem::mesh::mesh::print() const { +std::string specfem::mesh::mesh::print() const { - int n_elastic; - int n_acoustic; + std::ostringstream message; - Kokkos::parallel_reduce( - "specfem::mesh::mesh::print", specfem::kokkos::HostRange(0, this->nspec), - KOKKOS_CLASS_LAMBDA(const int ispec, int &n_elastic, int &n_acoustic) { - if (this->materials.material_index_mapping(ispec).type == - specfem::element::medium_tag::elastic) { - n_elastic++; - } else if (this->materials.material_index_mapping(ispec).type == - specfem::element::medium_tag::acoustic) { - n_acoustic++; - } - }, - n_elastic, n_acoustic); + int nspec = this->parameters.nspec; + int nglob = this->parameters.nglob; + int nspec_irregular = this->parameters.nspec_irregular; + int ngllx = this->parameters.ngllx; + int nglly = this->parameters.nglly; + int ngllz = this->parameters.ngllz; - std::ostringstream message; + int nacoustic = this->elements_types.nacoustic; + int nelastic = this->elements_types.nelastic; + int nporoelastic = this->elements_types.nporoelastic; - message - << "Spectral element information:\n" - << "------------------------------\n" - << "Total number of spectral elements : " << this->nspec << "\n" - << "Total number of spectral elements assigned to elastic material : " - << n_elastic << "\n" - << "Total number of spectral elements assigned to acoustic material : " - << n_acoustic << "\n" - << "Total number of geometric points : " << this->npgeo << "\n"; + // Print Mapping parameters + message << "3D Mesh information :\n" + << "------------------------------\n" + << "Total number of spectral elements: " << nspec << "\n" + << "Total number of global nodes: " << nglob << "\n" + << "Total number of irregular spectral elements: " << nspec_irregular + << "\n" + << "Total number of GLLX: " << ngllx << "\n" + << "Total number of GLLY: " << nglly << "\n" + << "Total number of GLLZ: " << ngllz << "\n" + << "\n" + << "Total number of acoustic spectral elements: " << nacoustic << "\n" + << "Total number of elastic spectral elements: " << nelastic << "\n" + << "Total number of poroelastic spectral elements: " << nporoelastic + << "\n" + << "------------------------------\n"; return message.str(); } diff --git a/src/mesh/dim3/parameters/parameters.cpp b/src/mesh/dim3/parameters/parameters.cpp index 24e1f460f..47e2e05ac 100644 --- a/src/mesh/dim3/parameters/parameters.cpp +++ b/src/mesh/dim3/parameters/parameters.cpp @@ -1,101 +1,107 @@ #include "mesh/dim3/parameters/parameters.hpp" #include +#include -void specfem::mesh::parameters::print() const { - std::cout << "Acoustic simulation:........................... " - << acoustic_simulation << std::endl; - std::cout << "Elastic simulation:............................ " - << elastic_simulation << std::endl; - std::cout << "Poroelastic simulation:........................ " - << poroelastic_simulation << std::endl; - std::cout << "Anisotropy:.................................... " << anisotropy - << std::endl; - std::cout << "Stacey ABC:.................................... " << stacey_abc - << std::endl; - std::cout << "PML ABC:....................................... " << pml_abc - << std::endl; - std::cout << "Approximate ocean load:........................ " - << approximate_ocean_load << std::endl; - std::cout << "Use mesh coloring:............................. " - << use_mesh_coloring << std::endl; - std::cout << "Number of dimensions:.......................... " << ndim - << std::endl; - std::cout << "Number of GLLX:................................ " << ngllx - << std::endl; - std::cout << "Number of GLLY:................................ " << nglly - << std::endl; - std::cout << "Number of GLLZ:................................ " << ngllz - << std::endl; - std::cout << "Number of GLL square:.......................... " << ngllsquare - << std::endl; - std::cout << "Number of processors:.......................... " << nproc - << std::endl; - std::cout << "Number of spectral elements:................... " << nspec - << std::endl; - std::cout << "Number of spectral elements for poroelastic:... " << nspec_poro - << std::endl; - std::cout << "Number of global nodes:........................ " << nglob - << std::endl; - std::cout << "Number of global nodes for ocean:.............. " << nglob_ocean - << std::endl; - std::cout << "Number of spectral elements for 2D bottom:..... " - << nspec2D_bottom << std::endl; - std::cout << "Number of spectral elements for 2D top:........ " << nspec2D_top - << std::endl; - std::cout << "Number of spectral elements for 2D xmin:....... " - << nspec2D_xmin << std::endl; - std::cout << "Number of spectral elements for 2D xmax:....... " - << nspec2D_xmax << std::endl; - std::cout << "Number of spectral elements for 2D ymin:....... " - << nspec2D_ymin << std::endl; - std::cout << "Number of spectral elements for 2D ymax:....... " - << nspec2D_ymax << std::endl; - std::cout << "Number of irregular spectral elements:......... " - << nspec_irregular << std::endl; - std::cout << "Number of neighbors:........................... " - << num_neighbors << std::endl; - std::cout << "Number of faces on the surface:................ " - << nfaces_surface << std::endl; - std::cout << "Number of absorbing boundary faces:............ " - << num_abs_boundary_faces << std::endl; - std::cout << "Number of free surface faces:.................. " - << num_free_surface_faces << std::endl; - std::cout << "Number of acoustic-elastic faces:.............. " - << num_coupling_ac_el_faces << std::endl; - std::cout << "Number of acoustic-poroelastic faces:.......... " - << num_coupling_ac_po_faces << std::endl; - std::cout << "Number of elastic-poroelastic faces:........... " - << num_coupling_el_po_faces << std::endl; - std::cout << "Number of poroelastic-elastic faces:........... " - << num_coupling_po_el_faces << std::endl; - std::cout << "Number of external mesh interfaces:............ " - << num_interfaces_ext_mesh << std::endl; - std::cout << "Maximum number of interfaces:.................. " - << max_nibool_interfaces_ext_mesh << std::endl; - std::cout << "Number of inner acoustic spectral elements:.... " - << nspec_inner_acoustic << std::endl; - std::cout << "Number of outer acoustic spectral elements:.... " - << nspec_outer_acoustic << std::endl; - std::cout << "Number of inner elastic spectral elements:..... " - << nspec_inner_elastic << std::endl; - std::cout << "Number of outer elastic spectral elements:..... " - << nspec_outer_elastic << std::endl; - std::cout << "Number of inner poroelastic spectral elements:. " - << nspec_inner_poroelastic << std::endl; - std::cout << "Number of outer poroelastic spectral elements:. " - << nspec_outer_poroelastic << std::endl; - std::cout << "Number of phase ispec acoustic: ............... " - << num_phase_ispec_acoustic << std::endl; - std::cout << "Number of phase ispec elastic:................. " - << num_phase_ispec_elastic << std::endl; - std::cout << "Number of phase ispec poroelastic:............. " - << num_phase_ispec_poroelastic << std::endl; - std::cout << "Number of colors inner acoustic:............... " - << num_colors_inner_acoustic << std::endl; - std::cout << "Number of colors outer acoustic:............... " - << num_colors_outer_acoustic << std::endl; - std::cout << "Number of colors inner elastic:................ " - << num_colors_inner_elastic << std::endl; - std::cout << "Number of colors outer elastic:................ " - << num_colors_outer_elastic << std::endl; +std::string +specfem::mesh::parameters::print() const { + + std::ostringstream message; + + message << "Acoustic simulation:........................... " + << acoustic_simulation << "\n"; + message << "Elastic simulation:............................ " + << elastic_simulation << "\n"; + message << "Poroelastic simulation:........................ " + << poroelastic_simulation << "\n"; + message << "Anisotropy:.................................... " << anisotropy + << "\n"; + message << "Stacey ABC:.................................... " << stacey_abc + << "\n"; + message << "PML ABC:....................................... " << pml_abc + << "\n"; + message << "Approximate ocean load:........................ " + << approximate_ocean_load << "\n"; + message << "Use mesh coloring:............................. " + << use_mesh_coloring << "\n"; + message << "Number of dimensions:.......................... " << ndim << "\n"; + message << "Number of GLLX:................................ " << ngllx + << "\n"; + message << "Number of GLLY:................................ " << nglly + << "\n"; + message << "Number of GLLZ:................................ " << ngllz + << "\n"; + message << "Number of GLL square:.......................... " << ngllsquare + << "\n"; + message << "Number of processors:.......................... " << nproc + << "\n"; + message << "Number of spectral elements:................... " << nspec + << "\n"; + message << "Number of spectral elements for poroelastic:... " << nspec_poro + << "\n"; + message << "Number of global nodes:........................ " << nglob + << "\n"; + message << "Number of global nodes for ocean:.............. " << nglob_ocean + << "\n"; + message << "Number of spectral elements for 2D bottom:..... " + << nspec2D_bottom << "\n"; + message << "Number of spectral elements for 2D top:........ " << nspec2D_top + << "\n"; + message << "Number of spectral elements for 2D xmin:....... " << nspec2D_xmin + << "\n"; + message << "Number of spectral elements for 2D xmax:....... " << nspec2D_xmax + << "\n"; + message << "Number of spectral elements for 2D ymin:....... " << nspec2D_ymin + << "\n"; + message << "Number of spectral elements for 2D ymax:....... " << nspec2D_ymax + << "\n"; + message << "Number of irregular spectral elements:......... " + << nspec_irregular << "\n"; + message << "Number of neighbors:........................... " << num_neighbors + << "\n"; + message << "Number of faces on the surface:................ " + << nfaces_surface << "\n"; + message << "Number of absorbing boundary faces:............ " + << num_abs_boundary_faces << "\n"; + message << "Number of free surface faces:.................. " + << num_free_surface_faces << "\n"; + message << "Number of acoustic-elastic faces:.............. " + << num_coupling_ac_el_faces << "\n"; + message << "Number of acoustic-poroelastic faces:.......... " + << num_coupling_ac_po_faces << "\n"; + message << "Number of elastic-poroelastic faces:........... " + << num_coupling_el_po_faces << "\n"; + message << "Number of poroelastic-elastic faces:........... " + << num_coupling_po_el_faces << "\n"; + message << "Number of external mesh interfaces:............ " + << num_interfaces_ext_mesh << "\n"; + message << "Maximum number of interfaces:.................. " + << max_nibool_interfaces_ext_mesh << "\n"; + message << "Number of inner acoustic spectral elements:.... " + << nspec_inner_acoustic << "\n"; + message << "Number of outer acoustic spectral elements:.... " + << nspec_outer_acoustic << "\n"; + message << "Number of inner elastic spectral elements:..... " + << nspec_inner_elastic << "\n"; + message << "Number of outer elastic spectral elements:..... " + << nspec_outer_elastic << "\n"; + message << "Number of inner poroelastic spectral elements:. " + << nspec_inner_poroelastic << "\n"; + message << "Number of outer poroelastic spectral elements:. " + << nspec_outer_poroelastic << "\n"; + message << "Number of phase ispec acoustic: ............... " + << num_phase_ispec_acoustic << "\n"; + message << "Number of phase ispec elastic:................. " + << num_phase_ispec_elastic << "\n"; + message << "Number of phase ispec poroelastic:............. " + << num_phase_ispec_poroelastic << "\n"; + message << "Number of colors inner acoustic:............... " + << num_colors_inner_acoustic << "\n"; + message << "Number of colors outer acoustic:............... " + << num_colors_outer_acoustic << "\n"; + message << "Number of colors inner elastic:................ " + << num_colors_inner_elastic << "\n"; + message << "Number of colors outer elastic:................ " + << num_colors_outer_elastic << "\n"; + + return message.str(); } diff --git a/src/mesh/dim3/partial_derivatives/partial_derivatives.cpp b/src/mesh/dim3/partial_derivatives/partial_derivatives.cpp index f9179f2fc..153076464 100644 --- a/src/mesh/dim3/partial_derivatives/partial_derivatives.cpp +++ b/src/mesh/dim3/partial_derivatives/partial_derivatives.cpp @@ -1,30 +1,39 @@ #include "mesh/dim3/partial_derivatives/partial_derivatives.hpp" #include +#include -void specfem::mesh::partial_derivatives::print() +std::string +specfem::mesh::partial_derivatives::print() const { - std::cout << "Partial parameters:\n" - << "------------------------------\n" - << "Number of spectral elements: " << nspec << "\n" - << "Number of GLLX:............. " << ngllx << "\n" - << "Number of GLLY:............. " << nglly << "\n" - << "Number of GLLZ:............. " << ngllz << "\n" - << "------------------------------\n"; + std::ostringstream message; + message << "Partial parameters:\n" + << "------------------------------\n" + << "Number of spectral elements: " << nspec << "\n" + << "Number of GLLX:............. " << ngllx << "\n" + << "Number of GLLY:............. " << nglly << "\n" + << "Number of GLLZ:............. " << ngllz << "\n" + << "------------------------------\n"; + return message.str(); } -void specfem::mesh::partial_derivatives::print( +std::string +specfem::mesh::partial_derivatives::print( int ispec, int igllx, int iglly, int igllz) const { - std::cout << "Partial parameters for spectral element " << ispec << ":\n" - << "------------------------------\n" - << "dxid (x,y,z): " << xix(ispec, igllx, iglly, igllz) << "\n" - << "xiy:..... " << xiy(ispec, igllx, iglly, igllz) << "\n" - << "xiz:..... " << xiz(ispec, igllx, iglly, igllz) << "\n" - << "etax:.... " << etax(ispec, igllx, iglly, igllz) << "\n" - << "etay:.... " << etay(ispec, igllx, iglly, igllz) << "\n" - << "etaz:.... " << etaz(ispec, igllx, iglly, igllz) << "\n" - << "gammax:.. " << gammax(ispec, igllx, iglly, igllz) << "\n" - << "gammay:.. " << gammay(ispec, igllx, iglly, igllz) << "\n" - << "gammaz:.. " << gammaz(ispec, igllx, iglly, igllz) << "\n" - << "jacobian: " << jacobian(ispec, igllx, iglly, igllz) << "\n" - << "------------------------------\n"; + + std::ostringstream message; + message << "Partial parameters for spectral element " << ispec << ":\n" + << "------------------------------\n" + << "dxid (x,y,z): " << xix(ispec, igllx, iglly, igllz) << "\n" + << "xiy:..... " << xiy(ispec, igllx, iglly, igllz) << "\n" + << "xiz:..... " << xiz(ispec, igllx, iglly, igllz) << "\n" + << "etax:.... " << etax(ispec, igllx, iglly, igllz) << "\n" + << "etay:.... " << etay(ispec, igllx, iglly, igllz) << "\n" + << "etaz:.... " << etaz(ispec, igllx, iglly, igllz) << "\n" + << "gammax:.. " << gammax(ispec, igllx, iglly, igllz) << "\n" + << "gammay:.. " << gammay(ispec, igllx, iglly, igllz) << "\n" + << "gammaz:.. " << gammaz(ispec, igllx, iglly, igllz) << "\n" + << "jacobian: " << jacobian(ispec, igllx, iglly, igllz) << "\n" + << "------------------------------\n"; + + return message.str(); } diff --git a/tests/unit-tests/mesh/materials/materials.cpp b/tests/unit-tests/mesh/materials/materials.cpp index ed72a1e51..3e8b2c40a 100644 --- a/tests/unit-tests/mesh/materials/materials.cpp +++ b/tests/unit-tests/mesh/materials/materials.cpp @@ -64,8 +64,9 @@ const static std::unordered_map 8100001799.8227, 8099996400.35, 8099996400.35, 0.0, 9999, 9999) }) } }; -void check_test(const specfem::mesh::materials &computed, - const MaterialVectorType &expected) { +void check_test( + const specfem::mesh::materials &computed, + const MaterialVectorType &expected) { if (computed.n_materials != expected.size()) { throw std::runtime_error("Size of materials is not the same");