Issue 1601 - Implement check mesh resolution #1630

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open

lsawade wants to merge 10 commits into devel from issue-1601-squash-merge

+2,224 −17

core/specfem/assembly/CMakeLists.txt

-Original file line number
+Diff line change
@@ -1,4 +1,4 @@
-    file(GLOB_RECURSE COMPUTE_SOURCE_FILES "./*.cpp" "./**/*.cpp")
+    file(GLOB_RECURSE COMPUTE_SOURCE_FILES "*.cpp")
     add_library(assembly ${COMPUTE_SOURCE_FILES})
     add_library(specfem::assembly ALIAS assembly)
@@ Expand Down @@

core/specfem/assembly/assembly/dim2/assembly.cpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -2,6 +2,7 @@ @@
     #include "enumerations/interface.hpp"
     #include "specfem/io.hpp"
     #include "specfem/mesh.hpp"
+    #include "specfem/assembly/info.hpp"
     specfem::assembly::assembly<specfem::element::dimension_tag::dim2>::assembly(
         const specfem::mesh::mesh<dimension_tag> &mesh,
@@ Expand Down Expand Up @@
                                          this->edge_types, this->mesh };
       this->fields = { this->mesh, this->element_types, simulation };
+      this->info = { this->mesh, this->properties, this->element_types };
       if (allocate_boundary_values)
         this->boundary_values = { max_timesteps, this->mesh, this->element_types,
                                   this->boundaries };
@@ Expand Down Expand Up @@
               << "Total number of geometric points : "
               << this->mesh.element_grid.ngllz << "\n"
               << "Total number of distinct quadrature points : " << this->mesh.nglob
-              << "\n";
+              << "\n"
+              << this->info.string() << "\n";
       int total_elements = 0;
@@ Expand Down @@

core/specfem/assembly/assembly/dim2/assembly.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -7,6 +7,7 @@ @@
     #include "specfem/assembly/conforming_interfaces.hpp"
     #include "specfem/assembly/edge_types.hpp"
     #include "specfem/assembly/element_types.hpp"
+    #include "specfem/assembly/info.hpp"
     #include "specfem/assembly/fields.hpp"
     #include "specfem/assembly/jacobian_matrix.hpp"
     #include "specfem/assembly/kernels.hpp"
@@ Expand Down Expand Up @@
                            ///< the
                            ///< boundaries
+      /**
+       * @brief Info
+       *
+       */
+      specfem::assembly::Info<dimension_tag> info;
       ///@}
       /**
@@ Expand Down Expand Up @@
        *
        */
       void check_jacobian_matrix() const;
     };
     } // namespace specfem::assembly

core/specfem/assembly/assembly/dim3/assembly.cpp

-Original file line number
+Diff line change
@@ Expand Up @@
       // Currently done in the mesher!
       this->check_jacobian_matrix();
+      this->info = { this->mesh, this->properties, this->element_types };
       return;
     }
@@ Expand All @@
               << "  Total number of spectral elements             : "
               << this->mesh.nspec << "\n"
               << "  Total number of quadrature points per element : "
-              << this->mesh.element_grid.ngllz << "\n";
-      // << "Total number of distinct quadrature points    : "
-      // << this->mesh.nglob << "\n";
+              << this->mesh.element_grid.ngllz << "\n"
+              << this->info.string() << "\n";
       int total_elements = 0;
@@ Expand Down @@

core/specfem/assembly/assembly/dim3/assembly.hpp

-Original file line number
+Diff line change
@@ Expand Up / @@ -6,6 +6,7 @@ @@
     #include "specfem/assembly/compute_source_array.hpp"
     #include "specfem/assembly/conforming_interfaces.hpp"
     #include "specfem/assembly/fields.hpp"
+    #include "specfem/assembly/info.hpp"
     #include "specfem/assembly/jacobian_matrix.hpp"
     #include "specfem/assembly/kernels.hpp"
     #include "specfem/assembly/mesh.hpp"
@@ Expand Down Expand Up @@
                            ///< the
                            ///< boundaries
+      specfem::assembly::Info<dimension_tag>
+          info; ///< Information about the mesh and
+                ///< simulation
       ///@}
       /**
@@ Expand Down @@

core/specfem/assembly/info.cpp

-Original file line number
+Diff line change
@@ -0,0 +1,37 @@
+    #include "specfem/assembly/info.hpp"
+    template <>
+    std::string specfem::assembly::Info<specfem::element::dimension_tag::dim2>::string() const {
+      std::ostringstream oss;
+      oss << "Mesh Information (2D):\n";
+      oss << " Domain X: ............... [" << domain_bounds.x().min << ", " << domain_bounds.x().max << "]\n";
+      oss << " Domain Z: ............... [" << domain_bounds.z().min << ", " << domain_bounds.z().max << "]\n";
+      oss << " VP: ..................... [" << vp.min << ", " << vp.max << "]\n";
+      oss << " VS: ..................... [" << vs.min << ", " << vs.max << "]\n";
+      oss << " V: .......................[" << v.min << ", " << v.max << "]\n";
+      oss << " Rho: .................... [" << rho.min << ", " << rho.max << "]\n";
+      oss << " Element Size: ........... [" << element_size.min << ", " << element_size.max << "]\n";
+      oss << " GLL Distance: ........... [" << gll_distance.min << ", " << gll_distance.max << "]\n";
+      oss << " Largest Minimum Period: . " << largest_minimum_period << "\n";
+      oss << " Suggested Time Step: .... " << suggested_time_step << "\n";
+      return oss.str();
+    }
+    template <>
+    std::string specfem::assembly::Info<specfem::element::dimension_tag::dim3>::string() const {
+      std::ostringstream oss;
+      oss << "Mesh Information (3D):\n";
+      oss << " Domain X: ............... [" << domain_bounds.x().min << ", " << domain_bounds.x().max << "]\n";
+      oss << " Domain Y: ............... [" << domain_bounds.y().min << ", " << domain_bounds.y().max << "]\n";
+      oss << " Domain Z: ............... [" << domain_bounds.z().min << ", " << domain_bounds.z().max << "]\n";
+      oss << " VP: ..................... [" << vp.min << ", " << vp.max << "]\n";
+      oss << " VS: ..................... [" << vs.min << ", " << vs.max << "]\n";
+      oss << " V: .......................[" << v.min << ", " << v.max << "]\n";
+      oss << " Rho: .................... [" << rho.min << ", " << rho.max << "]\n";
+      oss << " Element Size: ........... [" << element_size.min << ", " << element_size.max << "]\n";
+      oss << " GLL Distance: ........... [" << gll_distance.min << ", " << gll_distance.max << "]\n";
+      oss << " Largest Minimum Period: . " << largest_minimum_period << "\n";
+      oss << " Suggested Time Step: .... " << suggested_time_step << "\n";
+      return oss.str();
+    }

core/specfem/assembly/info.hpp

-Original file line number
+Diff line change
@@ -0,0 +1,70 @@
+    #pragma once
+    #include "specfem/assembly/element_types.hpp"
+    #include "specfem/assembly/info/impl/bounding_box.hpp"
+    #include "specfem/assembly/info/impl/bounds.hpp"
+    #include "specfem/assembly/mesh.hpp"
+    #include "specfem/assembly/properties.hpp"
+    #include "specfem/element.hpp"
+    #include "specfem_setup.hpp"
+    namespace specfem::assembly {
+    /**
+     * @brief Computes and stores mesh statistics and numerical stability
+     * parameters.
+     *
+     * Analyzes the assembled mesh to extract spatial bounds, material property
+     * ranges, and element geometry statistics. Estimates the minimum resolvable
+     * period and suggests a time step based on the CFL condition.
+     *
+     * The minimum period estimation follows Komatitsch et al. (2005):
+     * "average number of points per minimum wavelength in an element should be
+     * around 5."
+     *
+     * @tparam DimensionTag Spatial dimension (dim2 or dim3)
+     *
+     * @note The minimum period is an empirical estimate, not a sharp cutoff.
+     *       Synthetics become progressively less accurate for shorter periods.
+     */
+    template <specfem::element::dimension_tag DimensionTag> struct Info {
+      constexpr static auto dimension_tag = DimensionTag; ///< Dimension tag
+      Info() = default;
+      /**
+       * @brief Construct mesh info by analyzing mesh geometry and material
+       * properties.
+       *
+       * @param mesh Assembled mesh containing element geometry and GLL points
+       * @param properties Material properties at all mesh points
+       * @param element_types Element classification by medium and property type
+       */
+      Info(const specfem::assembly::mesh<dimension_tag> &mesh,
+           const specfem::assembly::properties<dimension_tag> &properties,
+           const specfem::assembly::element_types<dimension_tag> &element_types);
+      info::impl::BoundingBox<dimension_tag> domain_bounds; ///< Spatial extent of
+                                                            ///< the mesh domain
+      info::impl::Bounds element_size; ///< Element corner-to-corner distances
+      info::impl::Bounds gll_distance; ///< Distances between adjacent GLL points
+      info::impl::Bounds vp;           ///< P-wave velocity range
+      info::impl::Bounds vs;           ///< S-wave velocity range
+      info::impl::Bounds v;            ///< Combined wave velocity range
+      info::impl::Bounds rho;          ///< Density range
+      info::impl::Bounds vp_vs_ratio;  ///< Vp/Vs ratio range
+      type_real suggested_time_step;    ///< Time step satisfying CFL condition
+      type_real largest_minimum_period; ///< Maximum of minimum resolvable periods
+                                        ///< across elements
+      /**
+       * @brief Generate formatted string representation of mesh statistics.
+       * @return Multi-line string with labeled mesh properties
+       */
+      std::string string() const;
+    };
+    } // namespace specfem::assembly
+    #include "specfem/assembly/info.tpp"

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Issue 1601 - Implement check mesh resolution #1630

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Issue 1601 - Implement check mesh resolution #1630

Are you sure you want to change the base?

Issue 1601 - Implement check mesh resolution #1630

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!