fix boundary classifier

README.md

@@ -54,7 +54,7 @@ Both SEM and FD applications are compiled.
 
 Unit tests only
 ```sh
-ctest -LE benchmark
+ctest -LE "benchmark|validation"
 ```
 
 Benchmarks only, results will be stored in results generated in `build/Benchmarking` as a json file.

scripts/environments/env_Pangea4.sh

@@ -4,4 +4,5 @@
 module purge
 module load PrgEnv-gnu python/3.11.6 cmake/3.27.2 craype-x86-milan
 module unload cray-libsci cray-mpich
+module load gcc-native/13.2.0 # required for latest version of funtides
 module list

src/main/fe/src/sem_proxy.cc

@@ -91,26 +91,35 @@ SEMproxy::SEMproxy(const SemProxyOptions& opt)
             m_localParams.ex, m_localParams.lx, m_localParams.ey,
             m_localParams.ly, m_localParams.ez, m_localParams.lz,
             isModelOnNodes, isElastic, m_localParams.origin_x,
-            m_localParams.origin_y, m_localParams.origin_z);
-        m_mesh = builder.getModel();
+            m_localParams.origin_y, m_localParams.origin_z,
+            m_localParams.global_lx, m_localParams.global_ly,
+            m_localParams.global_lz, m_localParams.global_origin_x,
+            m_localParams.global_origin_y, m_localParams.global_origin_z);
+        m_mesh = builder.getModel(freeSurface_);
         break;
       }
       case 2: {
         model::CartesianStructBuilder<float, int, 2> builder(
             m_localParams.ex, m_localParams.lx, m_localParams.ey,
             m_localParams.ly, m_localParams.ez, m_localParams.lz,
             isModelOnNodes, isElastic, m_localParams.origin_x,
-            m_localParams.origin_y, m_localParams.origin_z);
-        m_mesh = builder.getModel();
+            m_localParams.origin_y, m_localParams.origin_z,
+            m_localParams.global_lx, m_localParams.global_ly,
+            m_localParams.global_lz, m_localParams.global_origin_x,
+            m_localParams.global_origin_y, m_localParams.global_origin_z);
+        m_mesh = builder.getModel(freeSurface_);
         break;
       }
       case 3: {
         model::CartesianStructBuilder<float, int, 3> builder(
             m_localParams.ex, m_localParams.lx, m_localParams.ey,
             m_localParams.ly, m_localParams.ez, m_localParams.lz,
             isModelOnNodes, isElastic, m_localParams.origin_x,
-            m_localParams.origin_y, m_localParams.origin_z);
-        m_mesh = builder.getModel();
+            m_localParams.origin_y, m_localParams.origin_z,
+            m_localParams.global_lx, m_localParams.global_ly,
+            m_localParams.global_lz, m_localParams.global_origin_x,
+            m_localParams.global_origin_y, m_localParams.global_origin_z);
+        m_mesh = builder.getModel(freeSurface_);
         break;
       }
       default:
@@ -122,7 +131,7 @@ SEMproxy::SEMproxy(const SemProxyOptions& opt)
   {
     // Pass local params to unstructured builder (handles origin internally)
     model::CartesianUnstructBuilder<float, int> builder(m_localParams);
-    m_mesh = builder.getModel();
+    m_mesh = builder.getModel(freeSurface_);
   }
   else
   {
@@ -232,10 +241,6 @@ void SEMproxy::run()
   const float taper_delta = 0.015;
 
   // Initialize Solver with Partition Info & Compute Local Mass
-
-  bool freeSurface = freeSurface_;
-  m_mesh->setFreeSurfaceEnabled(freeSurface);
-
   m_solver->computeFEInit(*m_mesh, sponge_size, surface_sponge, taper_delta);
 
   // Synchronize Mass Matrix (Critical for DD)

src/model/mesh/api/include/builder.h

@@ -18,8 +18,14 @@ class ModelBuilderBase
 
   static constexpr int MAX_ORDER = 5;
 
-  virtual std::shared_ptr<model::ModelApi<FloatType, ScalarType>> getModel()
-      const = 0;
+  /**
+   * @brief Get the model instance.
+   * @param free_surface_on_top Indicates if the free surface is on top, else we
+   * use damping on the top boundary.
+   * @return A shared pointer to the model instance.
+   */
+  virtual std::shared_ptr<model::ModelApi<FloatType, ScalarType>> getModel(
+      bool free_surface_on_top) const = 0;
 };
 }  // namespace model
 

src/model/mesh/api/include/model.h

@@ -119,13 +119,6 @@ class ModelApi
   PROXY_HOST_DEVICE
   virtual BoundaryFlag boundaryType(ScalarType n) const = 0;
 
-  /**
-   * @brief Initialize boundary flags based on node positions
-   * @param free_surface_on_top If true, mark top (Z+) as Surface, else as
-   * Damping
-   */
-  virtual void initializeBoundaryFlags(bool free_surface_on_top) = 0;
-
   /**
    * @brief Get P-wave velocity at a node
    * @param n Node index
@@ -350,21 +343,6 @@ class ModelApi
   PROXY_HOST_DEVICE
   virtual bool isFreeSurface(ScalarType n) const = 0;
 
-  /**
-   * @brief Enable/disable free surface on top boundary
-   * @param enable True to enable free surface BC, false for absorbing
-   * everywhere
-   */
-  virtual void setFreeSurfaceEnabled(bool enable) = 0;
-
-  /**
-   * @brief Initialize free surface detection
-   *
-   * Marks nodes on top boundary (Z+) as free surface.
-   * Called once during mesh setup.
-   */
-  virtual void initFreeSurface() = 0;
-
   /**
    * @brief Get total number of faces
    */

src/model/mesh/impl/builder/cartesian/include/cartesian_struct_boundary_classifier.h

@@ -0,0 +1,102 @@
+#ifndef FUNTIDES_MODEL_MESH_IMPL_BUILDER_CARTESIAN_INCLUDE_CARTESIAN_STRUCT_BOUNDARY_CLASSIFIER_H_
+#define FUNTIDES_MODEL_MESH_IMPL_BUILDER_CARTESIAN_INCLUDE_CARTESIAN_STRUCT_BOUNDARY_CLASSIFIER_H_
+
+#pragma once
+
+#include <data_type.h>
+#include <model.h>
+#include <model_struct.h>
+
+#include <cmath>
+
+namespace model
+{
+/**
+ * @brief Classifies boundary flags for nodes of a structured Cartesian mesh.
+ *
+ * Nodes are classified against the *global* domain bounds so that nodes on
+ * interior partition boundaries (e.g. MPI subdomain edges) are not marked as
+ * physical boundaries.
+ *
+ * Classification rules:
+ *  - Not on any global face              → InteriorNode
+ *  - On the z_max global face AND
+ *    free_surface_on_top == true         → Surface
+ *  - On any other global face            → Damping
+ *
+ * @tparam FloatType Floating-point type for coordinates and bounds
+ * @tparam ScalarType Integer type used to cast BoundaryFlag values
+ * @tparam Order Polynomial order of the spectral elements
+ */
+template <typename FloatType, typename ScalarType, int Order>
+class CartesianStructBoundaryClassifier
+{
+ public:
+  CartesianStructBoundaryClassifier(FloatType x_min, FloatType x_max,
+                                    FloatType y_min, FloatType y_max,
+                                    FloatType z_min, FloatType z_max,
+                                    bool free_surface_on_top)
+      : x_min_(x_min),
+        x_max_(x_max),
+        y_min_(y_min),
+        y_max_(y_max),
+        z_min_(z_min),
+        z_max_(z_max),
+        free_surface_on_top_(free_surface_on_top)
+  {
+  }
+
+  /**
+   * @brief Classify every node of @p model against the global domain bounds.
+   *
+   * Tolerance is derived from the minimum grid spacing of the model.
+   *
+   * @param model Structured mesh whose nodes are to be classified
+   * @return VECTOR_INT_VIEW of size model.getNumberOfNodes() with BoundaryFlag
+   *         values
+   */
+  VECTOR_INT_VIEW classify(
+      const model::ModelStruct<FloatType, ScalarType, Order>& model) const
+  {
+    const int n_node = model.getNumberOfNodes();
+    const FloatType tol = model.getMinSpacing() * static_cast<FloatType>(1e-4);
+
+    auto boundaries_t = allocateVector<VECTOR_INT_VIEW>(n_node, "boundaries_t");
+
+    for (int n = 0; n < n_node; ++n)
+    {
+      const FloatType x = model.nodeCoord(n, 0);
+      const FloatType y = model.nodeCoord(n, 1);
+      const FloatType z = model.nodeCoord(n, 2);
+
+      const bool at_xmin = (fabs(x - x_min_) < tol);
+      const bool at_xmax = (fabs(x - x_max_) < tol);
+      const bool at_ymin = (fabs(y - y_min_) < tol);
+      const bool at_ymax = (fabs(y - y_max_) < tol);
+      const bool at_zmin = (fabs(z - z_min_) < tol);
+      const bool at_zmax = (fabs(z - z_max_) < tol);
+
+      const bool on_boundary =
+          at_xmin || at_xmax || at_ymin || at_ymax || at_zmin || at_zmax;
+
+      if (!on_boundary)
+        boundaries_t(n) = static_cast<ScalarType>(BoundaryFlag::InteriorNode);
+      else if (free_surface_on_top_ && at_zmax)
+        boundaries_t(n) = static_cast<ScalarType>(BoundaryFlag::Surface);
+      else
+        boundaries_t(n) = static_cast<ScalarType>(BoundaryFlag::Damping);
+    }
+
+    return boundaries_t;
+  }
+
+ private:
+  FloatType x_min_, x_max_;
+  FloatType y_min_, y_max_;
+  FloatType z_min_, z_max_;
+  bool free_surface_on_top_;
+};
+
+}  // namespace model
+
+#endif  // FUNTIDES_MODEL_MESH_IMPL_BUILDER_CARTESIAN_INCLUDE_CARTESIAN_STRUCT_BOUNDARY_CLASSIFIER_H_

src/model/mesh/impl/builder/cartesian/include/cartesian_struct_builder.h

@@ -6,6 +6,8 @@
 #include <builder.h>
 #include <model_struct.h>
 
+#include "cartesian_struct_boundary_classifier.h"
+
 namespace model
 {
 template <typename FloatType, typename ScalarType, int Order>
@@ -16,11 +18,10 @@ class CartesianStructBuilder : public ModelBuilderBase<FloatType, ScalarType>
                          FloatType ly, ScalarType ez, FloatType lz,
                          bool isModelOnNodes, bool isElastic,
                          FloatType ox = 0.0, FloatType oy = 0.0,
-                         FloatType oz = 0.0,
-                         // AJOUTER : paramètres globaux (optionnels)
-                         FloatType global_lx = -1.0, FloatType global_ly = -1.0,
-                         FloatType global_lz = -1.0, FloatType global_ox = 0.0,
-                         FloatType global_oy = 0.0, FloatType global_oz = 0.0)
+                         FloatType oz = 0.0, FloatType global_lx = -1.0,
+                         FloatType global_ly = -1.0, FloatType global_lz = -1.0,
+                         FloatType global_ox = 0.0, FloatType global_oy = 0.0,
+                         FloatType global_oz = 0.0)
       : ex_(ex),
         ey_(ey),
         ez_(ez),
@@ -43,8 +44,8 @@ class CartesianStructBuilder : public ModelBuilderBase<FloatType, ScalarType>
 
   ~CartesianStructBuilder() = default;
 
-  std::shared_ptr<model::ModelApi<FloatType, ScalarType>> getModel()
-      const override
+  std::shared_ptr<model::ModelApi<FloatType, ScalarType>> getModel(
+      bool free_surface_on_top) const override
   {
     model::ModelStructData<FloatType, ScalarType> data;
     data.ex_ = ex_;
@@ -61,42 +62,17 @@ class CartesianStructBuilder : public ModelBuilderBase<FloatType, ScalarType>
 
     auto temp_model = model::ModelStruct<FloatType, ScalarType, Order>(data);
 
-    const int n_node = temp_model.getNumberOfNodes();
-    FloatType tol = temp_model.getMinSpacing() * 1e-4;
-
-    FloatType x_min = global_ox_, x_max = global_ox_ + global_lx_;
-    FloatType y_min = global_oy_, y_max = global_oy_ + global_ly_;
-    FloatType z_min = global_oz_, z_max = global_oz_ + global_lz_;
-
-    auto boundaries_t =
-        allocateVector<VECTOR_REAL_VIEW>(n_node, "boundaries_t");
-
-    for (int n = 0; n < n_node; ++n)
-    {
-      FloatType x = temp_model.nodeCoord(n, 0);
-      FloatType y = temp_model.nodeCoord(n, 1);
-      FloatType z = temp_model.nodeCoord(n, 2);
-
-      bool at_xmin = (fabs(x - x_min) < tol);
-      bool at_xmax = (fabs(x - x_max) < tol);
-      bool at_ymin = (fabs(y - y_min) < tol);
-      bool at_ymax = (fabs(y - y_max) < tol);
-      bool at_zmin = (fabs(z - z_min) < tol);
-      bool at_zmax = (fabs(z - z_max) < tol);
-
-      bool on_boundary =
-          at_xmin || at_xmax || at_ymin || at_ymax || at_zmin || at_zmax;
-
-      if (!on_boundary)
-        boundaries_t(n) = static_cast<FloatType>(BoundaryFlag::InteriorNode);
-      else if (at_zmax)
-        boundaries_t(n) = static_cast<FloatType>(BoundaryFlag::Surface);
-      else
-        boundaries_t(n) = static_cast<FloatType>(BoundaryFlag::Damping);
-    }
-
-    data.boundaries_t_ = boundaries_t;
-
+    data.boundaries_t_ =
+        CartesianStructBoundaryClassifier<FloatType, ScalarType, Order>(
+            global_ox_, global_ox_ + global_lx_, global_oy_,
+            global_oy_ + global_ly_, global_oz_, global_oz_ + global_lz_,
+            free_surface_on_top)
+            .classify(temp_model);
+
+    // -------------------------------------------------------------------------
+    // Construct model with local coordinates and dimensions, but use global
+    // boundaries for boundary classification.
+    // -------------------------------------------------------------------------
     auto model =
         std::make_shared<model::ModelStruct<FloatType, ScalarType, Order>>(
             data);