Adding Complex Boundaries

cli/new_command.py

@@ -24,7 +24,7 @@
 
 
 def ValidateSimName(sim_name):
-    pattern = re.compile("^[a-zA-Z]+[a-zA-Z0-9\-_]+$")
+    pattern = re.compile(r"^[a-zA-Z]+[a-zA-Z0-9\-_]+$")
     if not pattern.match(sim_name):
         Print.error('Error: simulation name "{}" is not valid.'.format(sim_name))
         Print.error("       Allowed characters are a-z A-Z 0-9 - and _")

cli/util.py

@@ -23,5 +23,5 @@
 def GetBinaryName():
     with open("CMakeLists.txt") as f:
         content = f.read()
-        return re.search("project\((.*)\)", content).group(1)
+        return re.search(r"project\((.*)\)", content).group(1)
 

cmake/external/ROOT.cmake

@@ -13,7 +13,11 @@ if(APPLE)
      "${DETECTED_OS_VERS}" MATCHES "^osx-11.7")
     execute_process(COMMAND bash "-c" "xcodebuild -version | sed -En 's/Xcode[[:space:]]+([0-9\.]*)/\\1/p'" OUTPUT_VARIABLE XCODE_VERS)
     message(STATUS "##### XCODE version: ${XCODE_VERS}")
-    if("${XCODE_VERS}" GREATER_EQUAL "15.0")
+    if("${XCODE_VERS}" GREATER_EQUAL "15.3")
+      message(STATUS "##### Using ROOT builds for XCODE 15.3")
+      set(ROOT_TAR_FILE root_v6.30.06_cxx17_python3.9_osx-xcode-15.3-${DETECTED_ARCH}.tar.gz)
+      set(ROOT_SHA_KEY osx-xcode-15.3-${DETECTED_ARCH}-ROOT)
+    elseif("${XCODE_VERS}" GREATER_EQUAL "15.0")
       message(STATUS "##### Using ROOT builds for XCODE 15.2")
       set(ROOT_TAR_FILE root_v6.30.02_cxx17_python3.9_osx-xcode-15.2-${DETECTED_ARCH}.tar.gz)
       set(ROOT_SHA_KEY osx-xcode-15.2-${DETECTED_ARCH}-ROOT)

cmake/external/SHA256Digests.cmake

@@ -58,6 +58,10 @@ SET(osx-14.2-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e8055
 SET(osx-14.2-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4)
 SET(osx-14.3-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e80554f7fa617ccb)
 SET(osx-14.3-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4)
+SET(osx-14.4-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e80554f7fa617ccb)
+SET(osx-14.4-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4)
+SET(osx-14.5-i386-ParaView d2e89df30ab0e2729b28539de37753e09c061b4c434b0a68e80554f7fa617ccb)
+SET(osx-14.5-arm64-ParaView 5e89b785ac0c56bbca31e2ae101a8986953e90227fe42a1ef7adcc70a0ff6fc4)
 SET(osx-xcode-13.1-i386-ROOT be97dd72022c8d082fbe4394f18b55c4920f20b138cfff1b5fc2b41d397ac203)
 SET(osx-xcode-13.1-arm64-ROOT 0a55b91c6df42d152b7943912e134f05c3872a73e73fcc129ee87fd847240ec8)
 SET(osx-xcode-14.1-i386-ROOT 001311608512b24535bb9710b8baf006bd00e9b0595fd6bdf900d28b1e22c395)
@@ -68,6 +72,8 @@ SET(osx-xcode-14.3-i386-ROOT c4f4c7a9848121ae19f3fe3c40c5812813ed42d9dfa7f1948f6
 SET(osx-xcode-14.3-arm64-ROOT 6ffcfaadd4de627cffde2e42fded9175b730769071396d7abdabeefc2714839a)
 SET(osx-xcode-15.2-i386-ROOT 3dea500451dc2c84d150a7d5f4b27d2cba3e867d7c775b744ac898f9f6e51249)
 SET(osx-xcode-15.2-arm64-ROOT 677429937f13a44a088b71851f0af111509bf81a9a25a299c905b6c9f7249587)
+SET(osx-xcode-15.3-i386-ROOT 52cef545acf679c3e491c9ef14e18edee454c90309b7fcf5f36fc2e659bdbde5)
+SET(osx-xcode-15.3-arm64-ROOT 3ee0885c329fab9a965c278cc491db82dcbde3999a96ac80d6f5145fb6794d1e)
 
 SET(ubuntu-18.04-Libroadrunner bf9293c1c95d0b65227bd7e08c0048116ba851bcec5028ef72ea13762ef79276)
 SET(ubuntu-18.04-ParaView e3fd74b13e9a4086988c5104c6b8d95c56365d25b491706a8e72018d0e5c76f1)

src/core/diffusion/diffusion_grid.cc

@@ -35,6 +35,8 @@ std::string BoundaryTypeToString(const BoundaryConditionType& type) {
       return "Dirichlet";
     case BoundaryConditionType::kPeriodic:
       return "Periodic";
+    case BoundaryConditionType::kComplex:
+      return "Complex";
     default:
       return "unknown";
   }
@@ -52,6 +54,8 @@ BoundaryConditionType StringToBoundaryType(const std::string& type) {
     return BoundaryConditionType::kDirichlet;
   } else if (type == "Periodic") {
     return BoundaryConditionType::kPeriodic;
+  } else if (type == "Complex") {
+    return BoundaryConditionType::kComplex;
   } else {
     Log::Fatal("StringToBoundaryType", "Unknown boundary type: ", type);
     return BoundaryConditionType::kNeumann;
@@ -64,7 +68,9 @@ DiffusionGrid::DiffusionGrid(int substance_id,
     : dc_({{1 - dc, dc / 6, dc / 6, dc / 6, dc / 6, dc / 6, dc / 6}}),
       mu_(mu),
       resolution_(resolution),
-      boundary_condition_(std::make_unique<ConstantBoundaryCondition>(0.0)) {
+      boundary_condition_(std::make_unique<ConstantBoundaryCondition>(0.0)),
+      boundary_condition_neumann_(std::make_unique<ConstantBoundaryCondition>(0.0)),
+      boundary_condition_dirichlet_(std::make_unique<ConstantBoundaryCondition>(0.0)) {
   // Compatibility with new abstract interface
   SetContinuumId(substance_id);
   SetContinuumName(substance_name);
@@ -81,6 +87,24 @@ void DiffusionGrid::Initialize() {
                "')");
   }
 
+  if (bc_type_ == BoundaryConditionType::kComplex) {
+    auto* param = Simulation::GetActive()->GetParam();
+    complex_bct_ = {param->diffusion_bc_x_min,
+                    param->diffusion_bc_x_max,
+                    param->diffusion_bc_y_min,
+                    param->diffusion_bc_y_max,
+                    param->diffusion_bc_z_min,
+                    param->diffusion_bc_z_max};
+    for (int i=0; i<6; i++) {
+      if (complex_bct_[i] != "Neumann" && complex_bct_[i] != "Dirichlet" && 
+          complex_bct_[i] != "Periodic") {
+        Log::Fatal("DiffusionGrid::Initialize",
+                   "Complex boundaries: unknown boundary condition ", complex_bct_[i],
+                   ". Select 'Neumann', 'Dirichlet' or 'Periodic' for each side of boundary.");
+      }
+    }
+  }
+
   // Get neighbor grid dimensions
   auto* env = Simulation::GetActive()->GetEnvironment();
   auto bounds = env->GetDimensionThresholds();
@@ -145,6 +169,8 @@ void DiffusionGrid::Diffuse(real_t dt) {
     DiffuseWithNeumann(dt);
   } else if (bc_type_ == BoundaryConditionType::kPeriodic) {
     DiffuseWithPeriodic(dt);
+  } else if (bc_type_ == BoundaryConditionType::kComplex) {
+    DiffuseWithComplex(dt);
   } else {
     Log::Error(
         "EulerGrid::Diffuse", "Boundary condition of type '",
@@ -533,6 +559,16 @@ void DiffusionGrid::SetBoundaryCondition(
   boundary_condition_ = std::move(bc);
 }
 
+void DiffusionGrid::SetBoundaryCondition_neumann(
+    std::unique_ptr<BoundaryCondition> bc) {
+  boundary_condition_neumann_ = std::move(bc);
+}
+
+void DiffusionGrid::SetBoundaryCondition_dirichlet(
+    std::unique_ptr<BoundaryCondition> bc) {
+  boundary_condition_dirichlet_ = std::move(bc);
+}
+
 BoundaryCondition* DiffusionGrid::GetBoundaryCondition() const {
   return boundary_condition_.get();
 }

src/core/diffusion/diffusion_grid.h

@@ -37,7 +37,8 @@ enum class BoundaryConditionType {
   kNeumann,
   kOpenBoundaries,
   kClosedBoundaries,
-  kPeriodic
+  kPeriodic,
+  kComplex
 };
 
 enum class InteractionMode { kAdditive = 0, kExponential = 1, kLogistic = 2 };
@@ -121,6 +122,11 @@ class DiffusionGrid : public ScalarField {
   /// Compute the diffusion of the substance with Periodic boundary conditions
   /// (u_{-1} = {u_N-1}, u_{N} = u_{0}) for a given time step dt.
   virtual void DiffuseWithPeriodic(real_t dt) = 0;
+  /// Compute the diffusion of the substance with Complex boundary conditions,
+  /// where the boundary type can be Dirichlet, Neumann or Periodic depending on
+  /// the side of the environment.
+  virtual void DiffuseWithComplex(real_t dt) = 0;
+
 
   /// Calculates the gradient for each box in the diffusion grid.
   /// The gradient is calculated in each direction (x, y, z) as following:
@@ -323,6 +329,10 @@ class DiffusionGrid : public ScalarField {
   /// `ConstantBoundaryCondition`
   void SetBoundaryCondition(std::unique_ptr<BoundaryCondition> bc);
 
+  void SetBoundaryCondition_neumann(std::unique_ptr<BoundaryCondition> bc);
+
+  void SetBoundaryCondition_dirichlet(std::unique_ptr<BoundaryCondition> bc);
+
   /// @brief  Returns the boundary condition. Does not transfer ownership.
   /// @return const Pointer to the boundary condition
   BoundaryCondition* GetBoundaryCondition() const;
@@ -415,6 +425,14 @@ class DiffusionGrid : public ScalarField {
   BoundaryConditionType bc_type_ = BoundaryConditionType::kDirichlet;
   /// Object that implements the boundary conditions.
   std::unique_ptr<BoundaryCondition> boundary_condition_ = nullptr;
+  /// Object that implements the boundary condition for Neumann type when
+  /// using complex boundaries.
+  std::unique_ptr<BoundaryCondition> boundary_condition_neumann_ = nullptr;
+  /// Object that implements the boundary condition for Dirichlet type when
+  /// using complex boundaries.
+  std::unique_ptr<BoundaryCondition> boundary_condition_dirichlet_ = nullptr;
+  /// Stores the boundary condition types when using complex boundaries.
+  std::array<std::string, 6> complex_bct_ = {{""}};
   /// Flag to indicate if the grid is initialized
   bool initialized_ = false;
   /// Flag to indicate if we want to print information about the grid after

src/core/diffusion/euler_grid.cc

@@ -366,4 +366,150 @@ void EulerGrid::DiffuseWithPeriodic(real_t dt) {
   c1_.swap(c2_);
 }
 
+
+
+
+void EulerGrid::DiffuseWithComplex(real_t dt) {
+
+  const auto nx = resolution_;
+  const auto ny = resolution_;
+  const auto nz = resolution_;
+  const auto num_boxes = nx * ny * nz;
+
+  const real_t d = 1 - dc_[0];
+  const real_t dx = box_length_;
+
+  const auto sim_time = GetSimulatedTime();
+
+  constexpr size_t YBF = 16;
+#pragma omp parallel for collapse(2)
+  for (size_t yy = 0; yy < ny; yy += YBF) {
+    for (size_t z = 0; z < nz; z++) {
+      size_t ymax = yy + YBF;
+      if (ymax >= ny) {
+        ymax = ny;
+      }
+      for (size_t y = yy; y < ymax; y++) {
+        size_t x{0};
+        size_t c{0};
+        size_t l{0};
+        size_t r{0};
+        size_t n{0};
+        size_t s{0};
+        size_t b{0};
+        size_t t{0};
+        c = x + y * nx + z * nx * ny;
+#pragma omp simd
+        for (x = 0; x < nx; x++) {
+          l = c - 1;
+          r = c + 1;
+          n = c - nx;
+          s = c + nx;
+          b = c - nx * ny;
+          t = c + nx * ny;
+
+          // Clamp to avoid out of bounds access. Clamped values are initialized
+          // to a wrong value but will be overwritten by the boundary condition
+          // evaluation. All other values are correct.
+          real_t left{c1_[std::clamp(l, size_t{0}, num_boxes - 1)]};
+          real_t right{c1_[std::clamp(r, size_t{0}, num_boxes - 1)]};
+          real_t bottom{c1_[std::clamp(b, size_t{0}, num_boxes - 1)]};
+          real_t top{c1_[std::clamp(t, size_t{0}, num_boxes - 1)]};
+          real_t north{c1_[std::clamp(n, size_t{0}, num_boxes - 1)]};
+          real_t south{c1_[std::clamp(s, size_t{0}, num_boxes - 1)]};
+
+          real_t CF_x{2.0};
+          real_t CF_y{2.0};
+          real_t CF_z{2.0};
+
+          if (x == 0 || x == (nx - 1) || y == 0 || y == (ny - 1) || z == 0 ||
+              z == (nz - 1)) {
+
+            real_t real_x = grid_dimensions_[0] + x * box_length_;
+            real_t real_y = grid_dimensions_[2] + y * box_length_;
+            real_t real_z = grid_dimensions_[4] + z * box_length_;
+
+            real_t boundary_value_neumann = boundary_condition_neumann_->Evaluate(real_x, real_y, real_z, sim_time);
+            real_t boundary_value_dirichlet = boundary_condition_dirichlet_->Evaluate(real_x, real_y, real_z, sim_time);
+
+            if (x == 0) {
+              if (complex_bct_[0] == "Dirichlet") {
+                left = boundary_value_dirichlet;
+              } else if (complex_bct_[0] == "Neumann") {
+                left = - boundary_value_neumann * box_length_;
+                CF_x -= 1.0;
+              } else if (complex_bct_[0] == "Periodic") {
+                l = nx-1 + y * nx + z * nx * ny;
+                left = c1_[l];
+              }
+            } else if (x == (nx - 1)) {
+              if (complex_bct_[1] == "Dirichlet") {
+                right = boundary_value_dirichlet;
+              } else if (complex_bct_[1] == "Neumann") {
+                right = - boundary_value_neumann * box_length_;
+                CF_x -= 1.0;
+              } else if (complex_bct_[1] == "Periodic") {
+                r = 0 + y * nx + z * nx * ny;
+                right = c1_[r];
+              }
+            }
+
+	    if (y == 0) {
+              if (complex_bct_[2] == "Dirichlet") {
+                north = boundary_value_dirichlet;
+              } else if (complex_bct_[2] == "Neumann") {
+                north = - boundary_value_neumann * box_length_;
+                CF_y -= 1.0;
+              } else if (complex_bct_[2] == "Periodic") {
+                n = x + (ny-1) * nx + z * nx * ny;
+                north = c1_[n];
+              }
+            } else if (y == (ny - 1)) {
+              if (complex_bct_[3] == "Dirichlet") {
+                south = boundary_value_dirichlet;
+              } else if (complex_bct_[3] == "Neumann") {
+                south = - boundary_value_neumann * box_length_;
+                CF_y -= 1.0;
+              } else if (complex_bct_[3] == "Periodic") {
+                s = x + 0 * nx + z * nx * ny;
+                south = c1_[s];
+              }
+            }
+
+            if (z == 0) {
+              if (complex_bct_[4] == "Dirichlet") {
+                bottom = boundary_value_dirichlet;
+              } else if (complex_bct_[4] == "Neumann") {
+                bottom = - boundary_value_neumann * box_length_;
+                CF_z -= 1.0;
+              } else if (complex_bct_[4] == "Periodic") {
+                b = x + y * nx + (nz-1) * nx * ny;
+                bottom = c1_[b];
+              }
+            } else if (z == (nz - 1)) {
+              if (complex_bct_[5] == "Dirichlet") {
+                top = boundary_value_dirichlet;
+              } else if (complex_bct_[5] == "Neumann") {
+                top = - boundary_value_neumann * box_length_;
+                CF_z -= 1.0;
+              } else if (complex_bct_[5] == "Periodic") {
+                t = x + y * nx + 0 * nx * ny;
+                top = c1_[t];
+              }
+            }
+          }
+
+	  c2_[c] = c1_[c] * (1 - (mu_ * dt)) + ( (d * dt / (dx*dx) ) * (left + right - CF_x*c1_[c]) )
+                                             + ( (d * dt / (dx*dx) ) * (south + north - CF_y*c1_[c]) )
+                                             + ( (d * dt / (dx*dx) ) * (top + bottom - CF_z*c1_[c]) );
+
+          ++c;
+        }
+      }  // tile ny
+    }    // tile nz
+  }      // block ny
+  c1_.swap(c2_);
+}
+
+
 }  // namespace bdm

src/core/diffusion/euler_grid.h

@@ -53,6 +53,7 @@ class EulerGrid : public DiffusionGrid {
   void DiffuseWithDirichlet(real_t dt) override;
   void DiffuseWithNeumann(real_t dt) override;
   void DiffuseWithPeriodic(real_t dt) override;
+  void DiffuseWithComplex(real_t dt) override;
 
  private:
   BDM_CLASS_DEF_OVERRIDE(EulerGrid, 1);

src/core/param/param.cc

@@ -239,6 +239,12 @@ void Param::AssignFromConfig(const std::shared_ptr<cpptoml::table>& config) {
   BDM_ASSIGN_CONFIG_VALUE(max_bound, "simulation.max_bound");
   BDM_ASSIGN_CONFIG_VALUE(diffusion_boundary_condition,
                           "simulation.diffusion_boundary_condition");
+  BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_x_min, "simulation.diffusion_bc_x_min");
+  BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_x_max, "simulation.diffusion_bc_x_max");
+  BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_y_min, "simulation.diffusion_bc_y_min");
+  BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_y_max, "simulation.diffusion_bc_y_max");
+  BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_z_min, "simulation.diffusion_bc_z_min");
+  BDM_ASSIGN_CONFIG_VALUE(diffusion_bc_z_max, "simulation.diffusion_bc_z_max");
   BDM_ASSIGN_CONFIG_VALUE(diffusion_method, "simulation.diffusion_method");
   BDM_ASSIGN_CONFIG_VALUE(calculate_gradients,
                           "simulation.calculate_gradients");

src/core/param/param.h

@@ -243,6 +243,22 @@ struct Param {
   ///     diffusion_boundary_condition = "Neumann"
   std::string diffusion_boundary_condition = "Neumann";
 
+  /// Define the boundary condition of each side of the simulation space
+  /// when using complex boundaries [Neumann, Dirichlet, Periodic].
+  /// Default value: `"Neumann"`\n TOML config file:
+  ///
+  ///     [simulation]
+  ///     diffusion_bc_x_min = "Neumann"
+  ///     diffusion_bc_x_max = "Neumann"
+  ///                ....
+  std::string diffusion_bc_x_min = "Neumann";
+  std::string diffusion_bc_x_max = "Neumann";
+  std::string diffusion_bc_y_min = "Neumann";
+  std::string diffusion_bc_y_max = "Neumann";
+  std::string diffusion_bc_z_min = "Neumann";
+  std::string diffusion_bc_z_max = "Neumann";
+
+
   /// A string for determining diffusion type within the simulation space.
   /// Currently, only the method "euler" implementing a FTCS scheme is
   /// supported. See for instance here:

test/unit/core/diffusion_init_test.h

@@ -32,6 +32,7 @@ class TestGrid : public DiffusionGrid {
   void DiffuseWithDirichlet(real_t dt) override { return; };
   void DiffuseWithNeumann(real_t dt) override { return; };
   void DiffuseWithPeriodic(real_t dt) override { return; };
+  void DiffuseWithComplex(real_t dt) override { return; };
 
   void Swap() { c1_.swap(c2_); }