From f7b1bbfddf4f73b778e0bbc4d7042e64fa53cd13 Mon Sep 17 00:00:00 2001
From: stanbot8 <stanbot8@gmail.com>
Date: Wed, 4 Mar 2026 12:10:50 +0000
Subject: [PATCH] Skip FTCS solver for zero-diffusion grids

When the diffusion coefficient is 0 but decay is nonzero, apply
decay as a simple element-wise multiply instead of running the
full FTCS stencil. The Laplacian term vanishes with D=0, so the
stencil computation is wasted work.
---
 src/core/diffusion/diffusion_grid.h | 16 +++++++++----
 src/core/diffusion/euler_grid.cc    | 35 ++++++++++++++++++++++++-----
 2 files changed, 42 insertions(+), 9 deletions(-)

diff --git a/src/core/diffusion/diffusion_grid.h b/src/core/diffusion/diffusion_grid.h
index 0b5e74c18..fd2569a43 100644
--- a/src/core/diffusion/diffusion_grid.h
+++ b/src/core/diffusion/diffusion_grid.h
@@ -170,7 +170,7 @@ class DiffusionGrid : public ScalarField {
   [[deprecated("Use GetValue instead")]] real_t GetConcentration(
       const Real3& position) const {
     return GetValue(position);
-  };
+  }
   /// @brief  Get the concentration at specified index
   /// @param idx Flat index of the grid
   /// @return c1_[idx]
@@ -337,9 +337,7 @@ class DiffusionGrid : public ScalarField {
   void PrintInfo(std::ostream& out = std::cout);
 
   /// Print the information after initialization
-  void PrintInfoWithInitialization() {
-    print_info_with_initialization_ = true;
-  };
+  void PrintInfoWithInitialization() { print_info_with_initialization_ = true; }
 
   /// Returns if the grid has been initialized
   bool IsInitialized() const { return initialized_; }
@@ -371,6 +369,16 @@ class DiffusionGrid : public ScalarField {
                    const ParallelResizeVector<Real3>& old_gradients,
                    size_t old_resolution);
 
+  /// Apply decay without diffusion (used when diffusion coefficient is 0).
+  void ApplyDecayOnly(real_t dt) {
+    const real_t decay = 1 - mu_ * dt;
+#pragma omp parallel for
+    for (size_t i = 0; i < total_num_boxes_; i++) {
+      c2_[i] = c1_[i] * decay;
+    }
+    c1_.swap(c2_);
+  }
+
   /// The side length of each box
   real_t box_length_ = 0;
   /// the volume of each box
diff --git a/src/core/diffusion/euler_grid.cc b/src/core/diffusion/euler_grid.cc
index f8ac1b722..31894ba33 100644
--- a/src/core/diffusion/euler_grid.cc
+++ b/src/core/diffusion/euler_grid.cc
@@ -19,12 +19,17 @@
 namespace bdm {
 
 void EulerGrid::DiffuseWithClosedEdge(real_t dt) {
+  const real_t d = 1 - dc_[0];
+  if (d == 0) {
+    ApplyDecayOnly(dt);
+    return;
+  }
+
   const auto nx = resolution_;
   const auto ny = resolution_;
   const auto nz = resolution_;
 
   const real_t ibl2 = 1 / (box_length_ * box_length_);
-  const real_t d = 1 - dc_[0];
 
   constexpr size_t YBF = 16;
 #pragma omp parallel for collapse(2)
@@ -67,12 +72,17 @@ void EulerGrid::DiffuseWithClosedEdge(real_t dt) {
 }
 
 void EulerGrid::DiffuseWithOpenEdge(real_t dt) {
+  const real_t d = 1 - dc_[0];
+  if (d == 0) {
+    ApplyDecayOnly(dt);
+    return;
+  }
+
   const auto nx = resolution_;
   const auto ny = resolution_;
   const auto nz = resolution_;
 
   const real_t ibl2 = 1 / (box_length_ * box_length_);
-  const real_t d = 1 - dc_[0];
   std::array<int, 4> l;
 
   constexpr size_t YBF = 16;
@@ -155,12 +165,17 @@ void EulerGrid::DiffuseWithOpenEdge(real_t dt) {
 }
 
 void EulerGrid::DiffuseWithDirichlet(real_t dt) {
+  const real_t d = 1 - dc_[0];
+  if (d == 0) {
+    ApplyDecayOnly(dt);
+    return;
+  }
+
   const auto nx = resolution_;
   const auto ny = resolution_;
   const auto nz = resolution_;
 
   const real_t ibl2 = 1 / (box_length_ * box_length_);
-  const real_t d = 1 - dc_[0];
 
   const auto sim_time = GetSimulatedTime();
 
@@ -211,13 +226,18 @@ void EulerGrid::DiffuseWithDirichlet(real_t dt) {
 }
 
 void EulerGrid::DiffuseWithNeumann(real_t dt) {
+  const real_t d = 1 - dc_[0];
+  if (d == 0) {
+    ApplyDecayOnly(dt);
+    return;
+  }
+
   const size_t nx = resolution_;
   const size_t ny = resolution_;
   const size_t nz = resolution_;
   const size_t num_boxes = nx * ny * nz;
 
   const real_t ibl2 = 1 / (box_length_ * box_length_);
-  const real_t d = 1 - dc_[0];
 
   const auto sim_time = GetSimulatedTime();
 
@@ -302,12 +322,17 @@ void EulerGrid::DiffuseWithNeumann(real_t dt) {
 }
 
 void EulerGrid::DiffuseWithPeriodic(real_t dt) {
+  const real_t d = 1 - dc_[0];
+  if (d == 0) {
+    ApplyDecayOnly(dt);
+    return;
+  }
+
   const size_t nx = resolution_;
   const size_t ny = resolution_;
   const size_t nz = resolution_;
 
   const real_t dx = box_length_;
-  const real_t d = 1 - dc_[0];
 
   constexpr size_t YBF = 16;
 #pragma omp parallel for collapse(2)