From f00c8156841be009dbd4b7fdd38d621705cdfc4e Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 17:20:00 -0700 Subject: [PATCH 1/2] Expose DiffusionGrid buffer access to subclasses via protected accessors Adds four protected accessors on DiffusionGrid: - GetConcentrationPtr / GetScratchPtr : mutable pointers to the live and scratch concentration buffers - SwapBuffers : swap the live and scratch buffers at the end of a step - GetRawDiffusionCoefficient : returns D directly (the stored dc_[0] is 1 - D for stencil convenience) Lets solver subclasses implement their own Diffuse* methods without requiring a friend declaration on DiffusionGrid. --- src/core/diffusion/diffusion_grid.h | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/core/diffusion/diffusion_grid.h b/src/core/diffusion/diffusion_grid.h index 0b5e74c18..50b875c80 100644 --- a/src/core/diffusion/diffusion_grid.h +++ b/src/core/diffusion/diffusion_grid.h @@ -348,6 +348,20 @@ class DiffusionGrid : public ScalarField { /// can be calculated on the fly. void TurnOffGradientCalculation() { precompute_gradients_ = false; } + protected: + /// Mutable access to the concentration buffer. Subclass solvers read from + /// this array in their stencil step. + real_t* GetConcentrationPtr() { return c1_.data(); } + /// Mutable access to the scratch buffer where the next-timestep values are + /// written before being swapped into the live buffer via SwapBuffers(). + real_t* GetScratchPtr() { return c2_.data(); } + /// Swap the live and scratch buffers. Call at the end of a step after + /// writing all next-timestep values into the scratch buffer. + void SwapBuffers() { c1_.swap(c2_); } + /// Raw diffusion coefficient D. The stored dc_[0] is 1 - D for stencil + /// convenience; this accessor returns D directly. + real_t GetRawDiffusionCoefficient() const { return 1 - dc_[0]; } + private: friend class EulerGrid; friend class EulerDepletionGrid; From 417067b971645f88906ccc971c8e93bda3f6a74e Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 17:45:00 -0700 Subject: [PATCH 2/2] Convert EulerGrid and EulerDepletionGrid to use protected accessors Each Diffuse* method hoists local aliases at the top and the loops use those aliases. Mechanical rewrite: c1_[i] -> c1[i], mu_ -> mu, c1_.swap(c2_) -> SwapBuffers(), boundary_condition_->Evaluate(...) -> GetBoundaryCondition()->Evaluate(...). Stencil and loop structure are unchanged. With no subclass reaching into private members, the EulerGrid and EulerDepletionGrid friend declarations are removed. Existing diffusion and euler tests still pass. --- src/core/diffusion/diffusion_grid.h | 2 - src/core/diffusion/euler_depletion_grid.cc | 13 +- src/core/diffusion/euler_grid.cc | 186 ++++++++++++--------- 3 files changed, 116 insertions(+), 85 deletions(-) diff --git a/src/core/diffusion/diffusion_grid.h b/src/core/diffusion/diffusion_grid.h index 50b875c80..8b22b7543 100644 --- a/src/core/diffusion/diffusion_grid.h +++ b/src/core/diffusion/diffusion_grid.h @@ -363,8 +363,6 @@ class DiffusionGrid : public ScalarField { real_t GetRawDiffusionCoefficient() const { return 1 - dc_[0]; } private: - friend class EulerGrid; - friend class EulerDepletionGrid; friend class TestGrid; // class used for testing (e.g. initialization) void ParametersCheck(real_t dt); diff --git a/src/core/diffusion/euler_depletion_grid.cc b/src/core/diffusion/euler_depletion_grid.cc index 6cd816f1b..cc9f1adac 100644 --- a/src/core/diffusion/euler_depletion_grid.cc +++ b/src/core/diffusion/euler_depletion_grid.cc @@ -19,6 +19,9 @@ namespace bdm { void EulerDepletionGrid::ApplyDepletion(real_t dt) { + const size_t nboxes = GetNumBoxes(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); auto* sim = Simulation::GetActive(); const auto* rm = sim->GetResourceManager(); @@ -27,7 +30,7 @@ void EulerDepletionGrid::ApplyDepletion(real_t dt) { // want to continue to use c1 for the next step. Thus, we swap pointers here // (and again after the depletion). This is necessary because ApplyDepletion // is called after the diffusion of the EulerGrid (swaps pointer at the end). - std::swap(c1_, c2_); + std::swap(c1, c2); for (size_t s = 0; s < binding_substances_.size(); s++) { if (binding_coefficients_[s] == 0.0) { @@ -50,14 +53,14 @@ void EulerDepletionGrid::ApplyDepletion(real_t dt) { auto* depleting_concentration = rm->GetDiffusionGrid(binding_substances_[s])->GetAllConcentrations(); #pragma omp parallel for simd - for (size_t c = 0; c < total_num_boxes_; c++) { - c2_[c] -= - c1_[c] * binding_coefficients_[s] * depleting_concentration[c] * dt; + for (size_t c = 0; c < nboxes; c++) { + c2[c] -= + c1[c] * binding_coefficients_[s] * depleting_concentration[c] * dt; } } } // See comment above. - std::swap(c1_, c2_); + std::swap(c1, c2); } void EulerDepletionGrid::DiffuseWithClosedEdge(real_t dt) { diff --git a/src/core/diffusion/euler_grid.cc b/src/core/diffusion/euler_grid.cc index f8ac1b722..192f25d6a 100644 --- a/src/core/diffusion/euler_grid.cc +++ b/src/core/diffusion/euler_grid.cc @@ -19,12 +19,18 @@ namespace bdm { void EulerGrid::DiffuseWithClosedEdge(real_t dt) { - 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 size_t res = GetResolution(); + const real_t bl = GetBoxLength(); + const real_t dc0 = GetDiffusionCoefficients()[0]; + const real_t mu = GetDecayConstant(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); + const auto nx = res; + const auto ny = res; + const auto nz = res; + + const real_t ibl2 = 1 / (bl * bl); + const real_t d = 1 - dc0; constexpr size_t YBF = 16; #pragma omp parallel for collapse(2) @@ -55,24 +61,30 @@ void EulerGrid::DiffuseWithClosedEdge(real_t dt) { b = c - nx * ny; t = c + nx * ny; - c2_[c] = c1_[c] * (1 - mu_ * dt) + - (d * dt * ibl2) * - (c1_[c - 1] - 2 * c1_[c] + c1_[c + 1] + c1_[s] - - 2 * c1_[c] + c1_[n] + c1_[b] - 2 * c1_[c] + c1_[t]); + c2[c] = + c1[c] * (1 - mu * dt) + + (d * dt * ibl2) * (c1[c - 1] - 2 * c1[c] + c1[c + 1] + c1[s] - + 2 * c1[c] + c1[n] + c1[b] - 2 * c1[c] + c1[t]); } } // tile ny } // tile nz } // block ny - c1_.swap(c2_); + SwapBuffers(); } void EulerGrid::DiffuseWithOpenEdge(real_t dt) { - 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 size_t res = GetResolution(); + const real_t bl = GetBoxLength(); + const real_t dc0 = GetDiffusionCoefficients()[0]; + const real_t mu = GetDecayConstant(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); + const auto nx = res; + const auto ny = res; + const auto nz = res; + + const real_t ibl2 = 1 / (bl * bl); + const real_t d = 1 - dc0; std::array l; constexpr size_t YBF = 16; @@ -122,10 +134,10 @@ void EulerGrid::DiffuseWithOpenEdge(real_t dt) { t = c + nx * ny; } - c2_[c] = c1_[c] * (1 - mu_ * dt) + - (d * dt * ibl2) * - (0 - 2 * c1_[c] + c1_[c + 1] + c1_[s] - 2 * c1_[c] + - c1_[n] + c1_[b] - 2 * c1_[c] + c1_[t]); + c2[c] = + c1[c] * (1 - mu * dt) + + (d * dt * ibl2) * (0 - 2 * c1[c] + c1[c + 1] + c1[s] - 2 * c1[c] + + c1[n] + c1[b] - 2 * c1[c] + c1[t]); #pragma omp simd for (x = 1; x < nx - 1; x++) { ++c; @@ -133,34 +145,40 @@ void EulerGrid::DiffuseWithOpenEdge(real_t dt) { ++s; ++b; ++t; - c2_[c] = - c1_[c] * (1 - mu_ * dt) + - (d * dt * ibl2) * (c1_[c - 1] - 2 * c1_[c] + c1_[c + 1] + - l[0] * c1_[s] - 2 * c1_[c] + l[1] * c1_[n] + - l[2] * c1_[b] - 2 * c1_[c] + l[3] * c1_[t]); + c2[c] = c1[c] * (1 - mu * dt) + + (d * dt * ibl2) * (c1[c - 1] - 2 * c1[c] + c1[c + 1] + + l[0] * c1[s] - 2 * c1[c] + l[1] * c1[n] + + l[2] * c1[b] - 2 * c1[c] + l[3] * c1[t]); } ++c; ++n; ++s; ++b; ++t; - c2_[c] = c1_[c] * (1 - mu_ * dt) + - (d * dt * ibl2) * - (c1_[c - 1] - 2 * c1_[c] + 0 + c1_[s] - 2 * c1_[c] + - c1_[n] + c1_[b] - 2 * c1_[c] + c1_[t]); + c2[c] = + c1[c] * (1 - mu * dt) + + (d * dt * ibl2) * (c1[c - 1] - 2 * c1[c] + 0 + c1[s] - 2 * c1[c] + + c1[n] + c1[b] - 2 * c1[c] + c1[t]); } // tile ny } // tile nz } // block ny - c1_.swap(c2_); + SwapBuffers(); } void EulerGrid::DiffuseWithDirichlet(real_t dt) { - 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 size_t res = GetResolution(); + const real_t bl = GetBoxLength(); + const auto& gd = GetDimensions(); + const real_t dc0 = GetDiffusionCoefficients()[0]; + const real_t mu = GetDecayConstant(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); + const auto nx = res; + const auto ny = res; + const auto nz = res; + + const real_t ibl2 = 1 / (bl * bl); + const real_t d = 1 - dc0; const auto sim_time = GetSimulatedTime(); @@ -185,11 +203,11 @@ void EulerGrid::DiffuseWithDirichlet(real_t dt) { if (x == 0 || x == (nx - 1) || y == 0 || y == (ny - 1) || z == 0 || z == (nz - 1)) { // For all boxes on the boundary, we simply evaluate the boundary - real_t real_x = grid_dimensions_[0] + x * box_length_; - real_t real_y = grid_dimensions_[0] + y * box_length_; - real_t real_z = grid_dimensions_[0] + z * box_length_; - c2_[c] = - boundary_condition_->Evaluate(real_x, real_y, real_z, sim_time); + real_t real_x = gd[0] + x * bl; + real_t real_y = gd[0] + y * bl; + real_t real_z = gd[0] + z * bl; + c2[c] = GetBoundaryCondition()->Evaluate(real_x, real_y, real_z, + sim_time); } else { // For inner boxes, we compute the regular stencil update n = c - nx; @@ -197,27 +215,34 @@ void EulerGrid::DiffuseWithDirichlet(real_t dt) { b = c - nx * ny; t = c + nx * ny; - c2_[c] = c1_[c] * (1 - mu_ * dt) + - (d * dt * ibl2) * - (c1_[c - 1] - 2 * c1_[c] + c1_[c + 1] + c1_[s] - - 2 * c1_[c] + c1_[n] + c1_[b] - 2 * c1_[c] + c1_[t]); + c2[c] = c1[c] * (1 - mu * dt) + + (d * dt * ibl2) * + (c1[c - 1] - 2 * c1[c] + c1[c + 1] + c1[s] - 2 * c1[c] + + c1[n] + c1[b] - 2 * c1[c] + c1[t]); } ++c; } } // tile ny } // tile nz } // block ny - c1_.swap(c2_); + SwapBuffers(); } void EulerGrid::DiffuseWithNeumann(real_t dt) { - const size_t nx = resolution_; - const size_t ny = resolution_; - const size_t nz = resolution_; + const size_t res = GetResolution(); + const real_t bl = GetBoxLength(); + const auto& gd = GetDimensions(); + const real_t dc0 = GetDiffusionCoefficients()[0]; + const real_t mu = GetDecayConstant(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); + const size_t nx = res; + const size_t ny = res; + const size_t nz = res; 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 real_t ibl2 = 1 / (bl * bl); + const real_t d = 1 - dc0; const auto sim_time = GetSimulatedTime(); @@ -247,22 +272,21 @@ void EulerGrid::DiffuseWithNeumann(real_t dt) { // 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(c - 1, size_t{0}, num_boxes - 1)]}; - real_t right{c1_[std::clamp(c + 1, 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 left{c1[std::clamp(c - 1, size_t{0}, num_boxes - 1)]}; + real_t right{c1[std::clamp(c + 1, 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 center_factor{6.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_[0] + y * box_length_; - real_t real_z = grid_dimensions_[0] + z * box_length_; - real_t boundary_value = - -box_length_ * - boundary_condition_->Evaluate(real_x, real_y, real_z, sim_time); + real_t real_x = gd[0] + x * bl; + real_t real_y = gd[0] + y * bl; + real_t real_z = gd[0] + z * bl; + real_t boundary_value = -bl * GetBoundaryCondition()->Evaluate( + real_x, real_y, real_z, sim_time); if (x == 0) { left = boundary_value; @@ -289,25 +313,31 @@ void EulerGrid::DiffuseWithNeumann(real_t dt) { } } - c2_[c] = c1_[c] * (1 - mu_ * dt) + - (d * dt * ibl2) * (left + right + south + north + top + - bottom - center_factor * c1_[c]); + c2[c] = c1[c] * (1 - mu * dt) + + (d * dt * ibl2) * (left + right + south + north + top + + bottom - center_factor * c1[c]); ++c; } } // tile ny } // tile nz } // block ny - c1_.swap(c2_); + SwapBuffers(); } void EulerGrid::DiffuseWithPeriodic(real_t dt) { - 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]; + const size_t res = GetResolution(); + const real_t bl = GetBoxLength(); + const real_t dc0 = GetDiffusionCoefficients()[0]; + const real_t mu = GetDecayConstant(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); + const size_t nx = res; + const size_t ny = res; + const size_t nz = res; + + const real_t dx = bl; + const real_t d = 1 - dc0; constexpr size_t YBF = 16; #pragma omp parallel for collapse(2) @@ -354,16 +384,16 @@ void EulerGrid::DiffuseWithPeriodic(real_t dt) { } // Stencil update - c2_[c] = c1_[c] * (1 - (mu_ * dt)) + - ((d * dt / (dx * dx)) * (c1_[l] + c1_[r] + c1_[n] + c1_[s] + - c1_[t] + c1_[b] - 6.0 * c1_[c])); + c2[c] = c1[c] * (1 - (mu * dt)) + + ((d * dt / (dx * dx)) * (c1[l] + c1[r] + c1[n] + c1[s] + + c1[t] + c1[b] - 6.0 * c1[c])); ++c; } } // tile ny } // tile nz } // block ny - c1_.swap(c2_); + SwapBuffers(); } } // namespace bdm