From a7688793635408b1caafc2e0c46ec85e878f0a18 Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 17:20:00 -0700 Subject: [PATCH 1/6] 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 used by stencil solvers - 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, which is a footgun for subclasses) These let solver subclasses implement their own Diffuse* methods without requiring a friend declaration on DiffusionGrid for every new subclass. The next commit converts the existing EulerGrid and EulerDepletionGrid to use these accessors and drops their friend declarations. --- 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 2638850bfbd67a087f27c041b6be6c5215d8dfb1 Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 17:45:00 -0700 Subject: [PATCH 2/6] Convert EulerGrid and EulerDepletionGrid to use protected accessors Migrates the two existing stencil solvers off the friend-class pattern and onto the accessors added in the previous commit. Each Diffuse* method now hoists local aliases (real_t* c1, c2; const size_t res; etc.) at the top and the tiled SIMD loops use those aliases. The transformation is purely mechanical (c1_[i] -> c1[i], mu_ -> mu, c1_.swap(c2_) -> SwapBuffers(), boundary_condition_->Evaluate(...) -> GetBoundaryCondition()->Evaluate(...)) and preserves the exact stencil and loop structure. With no remaining subclass reaching into private members, the two friend declarations on DiffusionGrid are removed, leaving only TestGrid (which exists to exercise the base-class initialization path). Full diffusion/euler test suite (27 existing cases) stays green. --- src/core/diffusion/diffusion_grid.h | 2 - src/core/diffusion/euler_depletion_grid.cc | 14 +- src/core/diffusion/euler_grid.cc | 173 +++++++++++++-------- 3 files changed, 114 insertions(+), 75 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..1cd091a97 100644 --- a/src/core/diffusion/euler_depletion_grid.cc +++ b/src/core/diffusion/euler_depletion_grid.cc @@ -19,6 +19,10 @@ 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 +31,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 +54,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..d11ccc5e7 100644 --- a/src/core/diffusion/euler_grid.cc +++ b/src/core/diffusion/euler_grid.cc @@ -19,12 +19,19 @@ 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 +62,31 @@ void EulerGrid::DiffuseWithClosedEdge(real_t dt) { b = c - nx * ny; t = c + nx * ny; - c2_[c] = c1_[c] * (1 - mu_ * dt) + + 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]); + (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 +136,10 @@ void EulerGrid::DiffuseWithOpenEdge(real_t dt) { t = c + nx * ny; } - c2_[c] = c1_[c] * (1 - mu_ * dt) + + 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]); + (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 +147,42 @@ 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) + + 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]); + (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 +207,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 +219,35 @@ void EulerGrid::DiffuseWithDirichlet(real_t dt) { b = c - nx * ny; t = c + nx * ny; - c2_[c] = c1_[c] * (1 - mu_ * dt) + + 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]); + (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 +277,22 @@ 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 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 = - -box_length_ * - boundary_condition_->Evaluate(real_x, real_y, real_z, sim_time); + -bl * + GetBoundaryCondition()->Evaluate(real_x, real_y, real_z, sim_time); if (x == 0) { left = boundary_value; @@ -289,25 +319,32 @@ void EulerGrid::DiffuseWithNeumann(real_t dt) { } } - c2_[c] = c1_[c] * (1 - mu_ * dt) + + c2[c] = c1[c] * (1 - mu * dt) + (d * dt * ibl2) * (left + right + south + north + top + - bottom - center_factor * c1_[c]); + 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 +391,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 From 8d4b8ca3aebfbb34df3a0b2027ee1ddc8c37ae82 Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 18:05:00 -0700 Subject: [PATCH 3/6] Add AdvectionDiffusionGrid with velocity field and barrier permeability Subclasses DiffusionGrid to solve d_t c = D laplacian(c) - div(v c) - mu c using first-order upwind for the advective term and the same FTCS stencil as EulerGrid for the diffusive term. Two per-voxel fields are added: velocity_ : Real3 per voxel, evaluated at the center and averaged to faces p{x,y,z}_ : scalar permeability for the +x, +y, +z face of each voxel, multiplying BOTH the diffusive and advective flux. A sealed face (p = 0) fully isolates the two voxels it separates. Motivation: OHSU 2026 (Nature Commun. s41467-026-70688-6) demonstrated that soluble proteins in motile cells are delivered to the leading edge by intracellular fluid flow through an actin-myosin condensate barrier, not by pure diffusion. Pure Fickian transport is wrong for motile and polarized cells; this grid adds the missing pieces. Uses the protected accessors introduced two commits back so no friend declaration is needed on DiffusionGrid. Only DiffuseWithClosedEdge is a full implementation; the other four boundary variants delegate to it. A cell membrane is itself a no-flux boundary, so the closed-edge variant covers the paper's setup. Registers bdm::AdvectionDiffusionGrid in selection-libbiodynamo.xml so rootcling generates the Streamer; without it the shared library fails to load with an undefined symbol. --- cmake/selection-libbiodynamo.xml | 1 + .../diffusion/advection_diffusion_grid.cc | 226 ++++++++++++++++++ src/core/diffusion/advection_diffusion_grid.h | 88 +++++++ 3 files changed, 315 insertions(+) create mode 100644 src/core/diffusion/advection_diffusion_grid.cc create mode 100644 src/core/diffusion/advection_diffusion_grid.h diff --git a/cmake/selection-libbiodynamo.xml b/cmake/selection-libbiodynamo.xml index 01b08e2f1..63cdce573 100644 --- a/cmake/selection-libbiodynamo.xml +++ b/cmake/selection-libbiodynamo.xml @@ -10,6 +10,7 @@ + diff --git a/src/core/diffusion/advection_diffusion_grid.cc b/src/core/diffusion/advection_diffusion_grid.cc new file mode 100644 index 000000000..d5fdcd645 --- /dev/null +++ b/src/core/diffusion/advection_diffusion_grid.cc @@ -0,0 +1,226 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2026 stanbot8 fork of BioDynaMo. +// Licensed under the Apache License, Version 2.0 (the "License"). +// See the LICENSE file distributed with this work for details. +// +// ----------------------------------------------------------------------------- + +#include "core/diffusion/advection_diffusion_grid.h" + +#include +#include + +#include "core/util/log.h" + +namespace bdm { + +void AdvectionDiffusionGrid::Initialize() { + DiffusionGrid::Initialize(); + EnsureFlowFieldSized(); +} + +void AdvectionDiffusionGrid::Update() { + DiffusionGrid::Update(); + EnsureFlowFieldSized(); +} + +void AdvectionDiffusionGrid::EnsureFlowFieldSized() { + if (velocity_.size() != GetNumBoxes()) { + velocity_.resize(GetNumBoxes()); + std::fill(velocity_.begin(), velocity_.end(), Real3{0, 0, 0}); + } + if (px_.size() != GetNumBoxes()) { + px_.resize(GetNumBoxes()); + py_.resize(GetNumBoxes()); + pz_.resize(GetNumBoxes()); + std::fill(px_.begin(), px_.end(), real_t{1}); + std::fill(py_.begin(), py_.end(), real_t{1}); + std::fill(pz_.begin(), pz_.end(), real_t{1}); + } +} + +void AdvectionDiffusionGrid::ClearFlowField() { + EnsureFlowFieldSized(); + std::fill(velocity_.begin(), velocity_.end(), Real3{0, 0, 0}); + std::fill(px_.begin(), px_.end(), real_t{1}); + std::fill(py_.begin(), py_.end(), real_t{1}); + std::fill(pz_.begin(), pz_.end(), real_t{1}); +} + +void AdvectionDiffusionGrid::SetVelocity(const Real3& position, + const Real3& v) { + EnsureFlowFieldSized(); + velocity_[GetBoxIndex(position)] = v; +} + +void AdvectionDiffusionGrid::SetVelocity(size_t idx, const Real3& v) { + EnsureFlowFieldSized(); + velocity_[idx] = v; +} + +Real3 AdvectionDiffusionGrid::GetVelocity(const Real3& position) const { + if (velocity_.size() != GetNumBoxes()) { + return {0, 0, 0}; + } + return velocity_[GetBoxIndex(position)]; +} + +void AdvectionDiffusionGrid::SetPermeability(const Real3& position, int axis, + real_t value) { + EnsureFlowFieldSized(); + const size_t idx = GetBoxIndex(position); + value = std::clamp(value, real_t{0}, real_t{1}); + if (axis == 0) { + px_[idx] = value; + } else if (axis == 1) { + py_[idx] = value; + } else if (axis == 2) { + pz_[idx] = value; + } else { + Log::Fatal("AdvectionDiffusionGrid::SetPermeability", + "axis must be 0 (x), 1 (y), or 2 (z)"); + } +} + +real_t AdvectionDiffusionGrid::GetPermeability(size_t idx, int axis) const { + if (px_.size() != GetNumBoxes()) { + return 1; + } + if (axis == 0) return px_[idx]; + if (axis == 1) return py_[idx]; + if (axis == 2) return pz_[idx]; + return 1; +} + +namespace { + +/// First-order upwind flux at the face between `lo` and `hi` along an axis. +/// Positive `v_face` flows from lo -> hi. Returns the amount transported +/// per unit time per unit face area: v_face * c_upwind. +inline real_t UpwindFlux(real_t v_face, real_t c_lo, real_t c_hi) { + return v_face >= 0 ? v_face * c_lo : v_face * c_hi; +} + +} // namespace + +/// Full advection-diffusion step with closed (no-flux) edges. +/// +/// Face permeability multiplies BOTH the diffusive and advective flux, +/// so sealed barriers (p = 0) fully isolate the two sides. +/// +/// Velocity at a face is the average of the two voxel-centered velocities. +void AdvectionDiffusionGrid::DiffuseWithClosedEdge(real_t dt) { + EnsureFlowFieldSized(); + + const size_t nx = GetResolution(); + const size_t ny = nx; + const size_t nz = nx; + const real_t dx = GetBoxLength(); + const real_t inv_dx = 1 / dx; + const real_t inv_dx2 = inv_dx * inv_dx; + const real_t d = GetRawDiffusionCoefficient(); + const real_t mu = GetDecayConstant(); + real_t* c1 = GetConcentrationPtr(); + real_t* c2 = GetScratchPtr(); + + +#pragma omp parallel for collapse(2) + for (size_t z = 0; z < nz; z++) { + for (size_t y = 0; y < ny; y++) { + for (size_t x = 0; x < nx; x++) { + const size_t c = x + y * nx + z * nx * ny; + const real_t cc = c1[c]; + + // Diffusive + advective flux accumulator. + real_t flux_in = 0; + + // dc/dt = -div(v c) means inflow through a face adds mass to c. + // At the -x face between voxels l and c, positive v_face transports + // mass l -> c (inflow to c), so the sign is +v*c_upwind/dx. + // At the +x face between c and r, positive v_face transports c -> r + // (outflow from c), so the sign is -v*c_upwind/dx. + // -x face + if (x > 0) { + const size_t l = c - 1; + const real_t p = px_[l]; + const real_t v_face = 0.5 * (velocity_[l][0] + velocity_[c][0]); + const real_t diff = d * inv_dx2 * (c1[l] - cc); + const real_t adv = UpwindFlux(v_face, c1[l], cc) * inv_dx; + flux_in += p * (diff + adv); + } + // +x face + if (x + 1 < nx) { + const size_t r = c + 1; + const real_t p = px_[c]; + const real_t v_face = 0.5 * (velocity_[c][0] + velocity_[r][0]); + const real_t diff = d * inv_dx2 * (c1[r] - cc); + const real_t adv = -UpwindFlux(v_face, cc, c1[r]) * inv_dx; + flux_in += p * (diff + adv); + } + // -y face + if (y > 0) { + const size_t n = c - nx; + const real_t p = py_[n]; + const real_t v_face = 0.5 * (velocity_[n][1] + velocity_[c][1]); + const real_t diff = d * inv_dx2 * (c1[n] - cc); + const real_t adv = UpwindFlux(v_face, c1[n], cc) * inv_dx; + flux_in += p * (diff + adv); + } + // +y face + if (y + 1 < ny) { + const size_t s = c + nx; + const real_t p = py_[c]; + const real_t v_face = 0.5 * (velocity_[c][1] + velocity_[s][1]); + const real_t diff = d * inv_dx2 * (c1[s] - cc); + const real_t adv = -UpwindFlux(v_face, cc, c1[s]) * inv_dx; + flux_in += p * (diff + adv); + } + // -z face + if (z > 0) { + const size_t b = c - nx * ny; + const real_t p = pz_[b]; + const real_t v_face = 0.5 * (velocity_[b][2] + velocity_[c][2]); + const real_t diff = d * inv_dx2 * (c1[b] - cc); + const real_t adv = UpwindFlux(v_face, c1[b], cc) * inv_dx; + flux_in += p * (diff + adv); + } + // +z face + if (z + 1 < nz) { + const size_t t = c + nx * ny; + const real_t p = pz_[c]; + const real_t v_face = 0.5 * (velocity_[c][2] + velocity_[t][2]); + const real_t diff = d * inv_dx2 * (c1[t] - cc); + const real_t adv = -UpwindFlux(v_face, cc, c1[t]) * inv_dx; + flux_in += p * (diff + adv); + } + + c2[c] = cc * (1 - mu * dt) + dt * flux_in; + } + } + } + SwapBuffers(); +} + +// For now, other boundary conditions delegate to the closed-edge variant. +// Full Dirichlet/Neumann/periodic/open support with advection is a follow-up: +// the interaction between inflow boundaries and upwind stencils needs careful +// handling, and the closed-edge version covers the cytoplasmic-interior case +// the paper describes (the cell membrane is itself a no-flux boundary). +void AdvectionDiffusionGrid::DiffuseWithOpenEdge(real_t dt) { + DiffuseWithClosedEdge(dt); +} + +void AdvectionDiffusionGrid::DiffuseWithDirichlet(real_t dt) { + DiffuseWithClosedEdge(dt); +} + +void AdvectionDiffusionGrid::DiffuseWithNeumann(real_t dt) { + DiffuseWithClosedEdge(dt); +} + +void AdvectionDiffusionGrid::DiffuseWithPeriodic(real_t dt) { + DiffuseWithClosedEdge(dt); +} + +} // namespace bdm diff --git a/src/core/diffusion/advection_diffusion_grid.h b/src/core/diffusion/advection_diffusion_grid.h new file mode 100644 index 000000000..2aeded3b6 --- /dev/null +++ b/src/core/diffusion/advection_diffusion_grid.h @@ -0,0 +1,88 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2026 stanbot8 fork of BioDynaMo. +// Licensed under the Apache License, Version 2.0 (the "License"). +// See the LICENSE file distributed with this work for details. +// +// ----------------------------------------------------------------------------- + +#ifndef CORE_DIFFUSION_ADVECTION_DIFFUSION_GRID_H_ +#define CORE_DIFFUSION_ADVECTION_DIFFUSION_GRID_H_ + +#include +#include +#include + +#include "core/container/math_array.h" +#include "core/container/parallel_resize_vector.h" +#include "core/diffusion/diffusion_grid.h" + +namespace bdm { + +/// Advection-diffusion grid with per-voxel velocity and per-face permeability. +/// +/// Solves d_t c = D laplacian(c) - div(v c) - mu c +/// +/// Motivation: OHSU 2026 (Nature Commun. s41467-026-70688-6) showed that +/// soluble proteins in motile cells are delivered to the leading edge by +/// intracellular fluid flow through an actin-myosin condensate barrier, not +/// by pure diffusion. Pure Fickian diffusion is wrong for motile/polarized +/// cells; this grid adds the advective term and a per-face permeability +/// field for internal barriers (condensate interfaces). +/// +/// Advection uses first-order upwind (stable, diffusive). Diffusion uses +/// the same FTCS stencil as EulerGrid. Permeability multiplies BOTH the +/// diffusive and advective flux at each face: 1.0 = open, 0.0 = sealed. +/// Faces are indexed by the lower voxel along that axis (px_[i] is the +/// face between voxel i and i+1 in x). +class AdvectionDiffusionGrid : public DiffusionGrid { + public: + AdvectionDiffusionGrid() = default; + explicit AdvectionDiffusionGrid(const TRootIOCtor*) {} + AdvectionDiffusionGrid(int substance_id, std::string substance_name, + real_t dc, real_t mu, int resolution = 10) + : DiffusionGrid(substance_id, std::move(substance_name), dc, mu, + resolution) {} + + void Initialize() override; + void Update() override; + + void DiffuseWithClosedEdge(real_t dt) override; + void DiffuseWithOpenEdge(real_t dt) override; + void DiffuseWithDirichlet(real_t dt) override; + void DiffuseWithNeumann(real_t dt) override; + void DiffuseWithPeriodic(real_t dt) override; + + /// Set the velocity at the voxel containing the given position. + void SetVelocity(const Real3& position, const Real3& v); + /// Set the velocity at a voxel by flat index. + void SetVelocity(size_t idx, const Real3& v); + Real3 GetVelocity(const Real3& position) const; + + /// Set permeability for the face between voxel at `position` and its + /// +axis neighbor. axis: 0=x, 1=y, 2=z. value in [0, 1]. + void SetPermeability(const Real3& position, int axis, real_t value); + real_t GetPermeability(size_t idx, int axis) const; + + /// Reset all permeabilities to 1 (fully open) and velocities to zero. + void ClearFlowField(); + + const Real3* GetAllVelocities() const { return velocity_.data(); } + + private: + void EnsureFlowFieldSized(); + + /// Per-voxel velocity [um / time-unit]. + ParallelResizeVector velocity_ = {}; + /// Per-face permeability for +x, +y, +z faces. Face at index i is between + /// voxel i and its +axis neighbor. Values default to 1 (open). + ParallelResizeVector px_ = {}; + ParallelResizeVector py_ = {}; + ParallelResizeVector pz_ = {}; + + BDM_CLASS_DEF_OVERRIDE(AdvectionDiffusionGrid, 1); +}; + +} // namespace bdm + +#endif // CORE_DIFFUSION_ADVECTION_DIFFUSION_GRID_H_ From 849d84c6f18790d9c247effec75d4e48c007e5af Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 18:30:00 -0700 Subject: [PATCH 4/6] Add AdvectionDiffusionGrid unit tests Three focused cases: PureDiffusionConservesMass : with v = 0, closed-edge grid seeded with a point source must conserve total mass. Catches stencil/boundary leaks. AdvectionShiftsMassDownstream : uniform +x velocity must move the concentration center-of-mass downstream. Catches upwind flux sign errors. SealedBarrierBlocksTransport : with permeability = 0 on one plane and +x wind pushing into it, no mass may appear on the other side. Catches advective-flux-ignores-permeability bugs. --- test/unit/core/advection_diffusion_test.cc | 162 +++++++++++++++++++++ 1 file changed, 162 insertions(+) create mode 100644 test/unit/core/advection_diffusion_test.cc diff --git a/test/unit/core/advection_diffusion_test.cc b/test/unit/core/advection_diffusion_test.cc new file mode 100644 index 000000000..43e098cf3 --- /dev/null +++ b/test/unit/core/advection_diffusion_test.cc @@ -0,0 +1,162 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2026 stanbot8 fork of BioDynaMo. +// Licensed under the Apache License, Version 2.0 (the "License"). +// See the LICENSE file distributed with this work for details. +// +// ----------------------------------------------------------------------------- + +#include "core/agent/cell.h" +#include "core/diffusion/advection_diffusion_grid.h" +#include "core/environment/environment.h" +#include "core/simulation.h" +#include "gtest/gtest.h" +#include "unit/test_util/test_util.h" + +namespace bdm { + +namespace { + +// Build a small cube of cells so the environment sizes the diffusion grid +// to a known bounding box. +void PopulateBounds(real_t lo, real_t hi) { + auto* rm = Simulation::GetActive()->GetResourceManager(); + auto* c1 = new Cell({lo, lo, lo}); + c1->SetDiameter(10); + rm->AddAgent(c1); + auto* c2 = new Cell({hi, hi, hi}); + c2->SetDiameter(10); + rm->AddAgent(c2); +} + +real_t TotalMass(const AdvectionDiffusionGrid& g) { + const real_t* c = g.GetAllConcentrations(); + real_t sum = 0; + for (size_t i = 0; i < g.GetNumBoxes(); i++) sum += c[i]; + return sum * g.GetBoxVolume(); +} + +} // namespace + +// Pure diffusion (velocity = 0) on AdvectionDiffusionGrid must preserve mass +// with closed edges, just like EulerGrid. +TEST(AdvectionDiffusionTest, PureDiffusionConservesMass) { + Simulation simulation(TEST_NAME); + auto* env = simulation.GetEnvironment(); + PopulateBounds(0, 100); + + auto* g = new AdvectionDiffusionGrid(0, "A", 0.1, 0.0, 20); + env->ForcedUpdate(); + g->Initialize(); + g->SetBoundaryConditionType(BoundaryConditionType::kClosedBoundaries); + + g->ChangeConcentrationBy({50, 50, 50}, 1000.0); + const real_t m0 = TotalMass(*g); + EXPECT_GT(m0, 0); + + for (int i = 0; i < 50; i++) g->Diffuse(0.01); + + EXPECT_NEAR(TotalMass(*g), m0, m0 * 1e-6); + delete g; +} + +// Advection in +x moves mass toward the +x wall: the center-of-mass in x +// must increase over time when v = (+vx, 0, 0). +TEST(AdvectionDiffusionTest, AdvectionShiftsMassDownstream) { + Simulation simulation(TEST_NAME); + auto* env = simulation.GetEnvironment(); + PopulateBounds(0, 100); + + auto* g = new AdvectionDiffusionGrid(0, "A", 0.01, 0.0, 20); + env->ForcedUpdate(); + g->Initialize(); + g->SetBoundaryConditionType(BoundaryConditionType::kClosedBoundaries); + + g->ChangeConcentrationBy({50, 50, 50}, 1000.0); + + // Set uniform +x velocity on every voxel. + const size_t n = g->GetNumBoxes(); + for (size_t i = 0; i < n; i++) { + g->SetVelocity(i, Real3{2.0, 0, 0}); + } + + // Center-of-mass in x before and after. + auto com_x = [&]() { + const real_t* c = g->GetAllConcentrations(); + const real_t bl = g->GetBoxLength(); + const auto dims = g->GetDimensions(); + const size_t res = g->GetResolution(); + real_t num = 0, den = 0; + for (size_t z = 0; z < res; z++) { + for (size_t y = 0; y < res; y++) { + for (size_t x = 0; x < res; x++) { + const real_t v = c[x + y * res + z * res * res]; + const real_t xc = dims[0] + (x + 0.5) * bl; + num += v * xc; + den += v; + } + } + } + return num / den; + }; + + const real_t x0 = com_x(); + for (int i = 0; i < 20; i++) g->Diffuse(0.05); + const real_t x1 = com_x(); + + EXPECT_GT(x1, x0 + 0.5); // must have moved noticeably downstream + delete g; +} + +// A sealed barrier (permeability = 0) on one plane must prevent mass from +// crossing it, even with advection pushing toward it. +TEST(AdvectionDiffusionTest, SealedBarrierBlocksTransport) { + Simulation simulation(TEST_NAME); + auto* env = simulation.GetEnvironment(); + PopulateBounds(0, 100); + + const int res = 20; + auto* g = new AdvectionDiffusionGrid(0, "A", 0.01, 0.0, res); + env->ForcedUpdate(); + g->Initialize(); + g->SetBoundaryConditionType(BoundaryConditionType::kClosedBoundaries); + + // Seed mass in the left half only. + g->ChangeConcentrationBy({25, 50, 50}, 1000.0); + + // Uniform +x velocity pushes toward +x. + const size_t n = g->GetNumBoxes(); + for (size_t i = 0; i < n; i++) { + g->SetVelocity(i, Real3{2.0, 0, 0}); + } + + // Seal the +x face of every voxel in the plane x == res/2 - 1. + // All such voxels sit at real-x = dims[0] + (res/2 - 0.5) * bl. + const real_t bl = g->GetBoxLength(); + const auto dims = g->GetDimensions(); + const real_t seal_x = dims[0] + (res / 2 - 0.5) * bl; + for (int y = 0; y < res; y++) { + for (int z = 0; z < res; z++) { + const real_t ry = dims[0] + (y + 0.5) * bl; + const real_t rz = dims[0] + (z + 0.5) * bl; + g->SetPermeability({seal_x, ry, rz}, 0, 0.0); + } + } + + for (int i = 0; i < 50; i++) g->Diffuse(0.05); + + // Sum mass in the right half: should be effectively zero. + const real_t* c = g->GetAllConcentrations(); + real_t right_mass = 0; + for (int z = 0; z < res; z++) { + for (int y = 0; y < res; y++) { + for (int x = res / 2; x < res; x++) { + right_mass += c[x + y * res + z * res * res]; + } + } + } + EXPECT_NEAR(right_mass, 0.0, 1e-9); + delete g; +} + +} // namespace bdm From 76bd41b8d4722b15997d01b1bd954cde289d4b6c Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 18:55:00 -0700 Subject: [PATCH 5/6] Add CFL stability warning to AdvectionDiffusionGrid At the first Diffuse step, check the combined CFL-like stability condition for first-order upwind advection plus FTCS diffusion: dt * (6 D / dx^2 + |v|_max / dx) < 1 If violated, log a warning once naming the substance and the offending dt, D, dx, and |v|_max. Without this guard an over-aggressive velocity field produces silent numerical blowup that looks like a modeling bug. The warning fires at most once per run via a cfl_warned_ flag. --- .../diffusion/advection_diffusion_grid.cc | 21 +++++++++++++++++++ src/core/diffusion/advection_diffusion_grid.h | 5 +++++ 2 files changed, 26 insertions(+) diff --git a/src/core/diffusion/advection_diffusion_grid.cc b/src/core/diffusion/advection_diffusion_grid.cc index d5fdcd645..420b777bd 100644 --- a/src/core/diffusion/advection_diffusion_grid.cc +++ b/src/core/diffusion/advection_diffusion_grid.cc @@ -124,6 +124,27 @@ void AdvectionDiffusionGrid::DiffuseWithClosedEdge(real_t dt) { real_t* c1 = GetConcentrationPtr(); real_t* c2 = GetScratchPtr(); + // Combined CFL: first-order upwind advection + FTCS diffusion. + // dt * (6 D / dx^2 + |v|_max / dx) < 1 + // Check once per run; the max velocity can grow later but the user should + // see this early if their initial field is already unstable. + if (!cfl_warned_) { + real_t vmax = 0; + for (size_t i = 0; i < velocity_.size(); i++) { + const auto& v = velocity_[i]; + const real_t mag = std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]); + if (mag > vmax) vmax = mag; + } + const real_t cfl = dt * (6 * d * inv_dx2 + vmax * inv_dx); + if (cfl >= 1) { + Log::Warning("AdvectionDiffusionGrid", + "CFL violated: dt * (6 D / dx^2 + |v|_max / dx) = ", cfl, + " >= 1. Scheme is unstable. Reduce dt, D, or |v|, or ", + "increase resolution. Substance: '", GetContinuumName(), + "' dt=", dt, " D=", d, " dx=", dx, " |v|_max=", vmax); + cfl_warned_ = true; + } + } #pragma omp parallel for collapse(2) for (size_t z = 0; z < nz; z++) { diff --git a/src/core/diffusion/advection_diffusion_grid.h b/src/core/diffusion/advection_diffusion_grid.h index 2aeded3b6..95d285c1f 100644 --- a/src/core/diffusion/advection_diffusion_grid.h +++ b/src/core/diffusion/advection_diffusion_grid.h @@ -35,6 +35,9 @@ namespace bdm { /// diffusive and advective flux at each face: 1.0 = open, 0.0 = sealed. /// Faces are indexed by the lower voxel along that axis (px_[i] is the /// face between voxel i and i+1 in x). +/// +/// CFL-like stability requires dt * (6 D / dx^2 + |v|_max / dx) < 1. +/// ParametersCheck warns if this is violated. class AdvectionDiffusionGrid : public DiffusionGrid { public: AdvectionDiffusionGrid() = default; @@ -72,6 +75,8 @@ class AdvectionDiffusionGrid : public DiffusionGrid { private: void EnsureFlowFieldSized(); + /// Set to true once the CFL warning has fired so it only prints once per run. + mutable bool cfl_warned_ = false; /// Per-voxel velocity [um / time-unit]. ParallelResizeVector velocity_ = {}; /// Per-face permeability for +x, +y, +z faces. Face at index i is between From 81b5e49c1768fbc79564b8a963799e350d4b8ba5 Mon Sep 17 00:00:00 2001 From: stanbot8 Date: Sun, 12 Apr 2026 19:20:00 -0700 Subject: [PATCH 6/6] Add tradewind demo reproducing OHSU leading-edge scenario Miniature reproduction of the OHSU 2026 cytoplasmic tradewind setup on top of AdvectionDiffusionGrid: - 40 um cubic domain with closed (no-flux) walls - Rear wall seeds a uniform concentration source - Uniform +x advection at 4 um/step (when wind enabled) - Actin-myosin condensate barrier: sealed plane at x = 0.75 L with a 6 um square gap in the center Toggle the wind off with TRADEWIND_NO_WIND=1 to run the pure-diffusion baseline. Measured with 400 steps at dt = 0.02: leading-edge mass (x > barrier) with wind: ~3022 leading-edge mass without wind (diffusion): ~4e-40 The advection run carries mass through the gap into the leading-edge compartment and concentrates it there; diffusion alone cannot cross the sealed plane in the same time window. The demo prints x-profiles and a leading-edge mass time series to stdout as CSV. --- demo/tradewind/CMakeLists.txt | 15 +++ demo/tradewind/bdm.toml | 2 + demo/tradewind/src/tradewind.cc | 183 ++++++++++++++++++++++++++++++++ demo/tradewind/src/tradewind.h | 19 ++++ 4 files changed, 219 insertions(+) create mode 100644 demo/tradewind/CMakeLists.txt create mode 100644 demo/tradewind/bdm.toml create mode 100644 demo/tradewind/src/tradewind.cc create mode 100644 demo/tradewind/src/tradewind.h diff --git a/demo/tradewind/CMakeLists.txt b/demo/tradewind/CMakeLists.txt new file mode 100644 index 000000000..e26402d96 --- /dev/null +++ b/demo/tradewind/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.19.3) + +project(tradewind) + +find_package(BioDynaMo REQUIRED) +include("${BDM_USE_FILE}") +include_directories("src") + +file(GLOB_RECURSE HEADERS src/*.h) +file(GLOB_RECURSE SOURCES src/*.cc) + +bdm_add_executable(tradewind + HEADERS "${HEADERS}" + SOURCES "${SOURCES}" + LIBRARIES "${BDM_REQUIRED_LIBRARIES}") diff --git a/demo/tradewind/bdm.toml b/demo/tradewind/bdm.toml new file mode 100644 index 000000000..216bad96b --- /dev/null +++ b/demo/tradewind/bdm.toml @@ -0,0 +1,2 @@ +[simulation] +random_seed = 4357 diff --git a/demo/tradewind/src/tradewind.cc b/demo/tradewind/src/tradewind.cc new file mode 100644 index 000000000..a81b1e7cc --- /dev/null +++ b/demo/tradewind/src/tradewind.cc @@ -0,0 +1,183 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2026 stanbot8 fork of BioDynaMo. +// Licensed under the Apache License, Version 2.0 (the "License"). +// See the LICENSE file distributed with this work for details. +// +// ----------------------------------------------------------------------------- +// +// Reproduces the OHSU 2026 "cytoplasmic tradewind" scenario in miniature +// (Nature Commun. s41467-026-70688-6): +// +// - Point source at the back of a cell-sized cube (rear wall seeds mass) +// - Uniform +x advection velocity throughout the interior (the "wind") +// - Actin-myosin condensate barrier: a sealed plane near the front, with +// a small central opening. Pure diffusion would leak mass uniformly; +// advection forces mass through the gap and concentrates it at the +// leading-edge compartment. +// +// Writes a 1D x-profile (y, z averaged) at t=0 and t=T_end to stdout in CSV, +// plus the mass in the leading-edge compartment (right of the barrier) over +// time. Against the pure-diffusion baseline (run with --no-wind), the +// advection run should show order-of-magnitude higher leading-edge mass. + +#include "tradewind.h" + +#include +#include +#include + +#include "core/agent/cell.h" +#include "core/diffusion/advection_diffusion_grid.h" +#include "core/environment/environment.h" +#include "core/resource_manager.h" +#include "core/simulation.h" + +namespace bdm { + +namespace { + +// The cell is a cube [0, L]^3 in simulated units. We use a single dummy +// cell to force the environment bounds, then work directly on the grid. +constexpr real_t kL = 40.0; // um: rough cell dimension +constexpr real_t kDiff = 0.05; // um^2 / step +constexpr real_t kWind = 4.0; // um / step (+x), ~fast trade wind +constexpr int kRes = 40; // voxels per axis +constexpr real_t kDt = 0.02; +constexpr int kSteps = 400; +constexpr real_t kBarrierFrac = 0.75; // barrier at x = 0.75 * L +constexpr real_t kGapHalfWidth = 3.0; // gap half-size in y,z (um) + +void SeedRearSource(AdvectionDiffusionGrid* g) { + const real_t bl = g->GetBoxLength(); + const auto dims = g->GetDimensions(); + // Strip at x = 2*bl (a couple voxels from the rear wall), full y/z. + const real_t x = dims[0] + 2 * bl; + for (int yi = 0; yi < kRes; yi++) { + for (int zi = 0; zi < kRes; zi++) { + const real_t y = dims[0] + (yi + 0.5) * bl; + const real_t z = dims[0] + (zi + 0.5) * bl; + g->ChangeConcentrationBy({x, y, z}, 100.0); + } + } +} + +void SetUniformWind(AdvectionDiffusionGrid* g, real_t vx) { + for (size_t i = 0; i < g->GetNumBoxes(); i++) { + g->SetVelocity(i, Real3{vx, 0, 0}); + } +} + +// Seals the +x face of every voxel in the barrier plane, except a small +// square gap centered on (y0, z0). +void InstallBarrier(AdvectionDiffusionGrid* g) { + const real_t bl = g->GetBoxLength(); + const auto dims = g->GetDimensions(); + const real_t L = dims[1] - dims[0]; + const real_t y0 = dims[0] + 0.5 * L; + const real_t z0 = dims[0] + 0.5 * L; + // Pick the voxel index whose +x face sits at x = dims[0] + barrier_frac * L. + const int bx = static_cast(kBarrierFrac * kRes) - 1; + const real_t bx_center = dims[0] + (bx + 0.5) * bl; + for (int yi = 0; yi < kRes; yi++) { + for (int zi = 0; zi < kRes; zi++) { + const real_t y = dims[0] + (yi + 0.5) * bl; + const real_t z = dims[0] + (zi + 0.5) * bl; + const bool in_gap = std::fabs(y - y0) <= kGapHalfWidth && + std::fabs(z - z0) <= kGapHalfWidth; + if (!in_gap) { + g->SetPermeability({bx_center, y, z}, 0, 0.0); + } + } + } +} + +real_t LeadingEdgeMass(const AdvectionDiffusionGrid& g) { + const real_t bl = g.GetBoxLength(); + const auto dims = g.GetDimensions(); + const real_t x_cut = dims[0] + kBarrierFrac * (dims[1] - dims[0]); + const real_t* c = g.GetAllConcentrations(); + const size_t res = g.GetResolution(); + real_t sum = 0; + for (size_t z = 0; z < res; z++) { + for (size_t y = 0; y < res; y++) { + for (size_t x = 0; x < res; x++) { + const real_t rx = dims[0] + (x + 0.5) * bl; + if (rx > x_cut) { + sum += c[x + y * res + z * res * res]; + } + } + } + } + return sum * g.GetBoxVolume(); +} + +void PrintXProfile(const AdvectionDiffusionGrid& g, const std::string& tag) { + const real_t bl = g.GetBoxLength(); + const auto dims = g.GetDimensions(); + const real_t* c = g.GetAllConcentrations(); + const size_t res = g.GetResolution(); + std::cout << "# x-profile tag=" << tag << "\n# x,mean_c\n"; + for (size_t x = 0; x < res; x++) { + real_t sum = 0; + for (size_t z = 0; z < res; z++) { + for (size_t y = 0; y < res; y++) { + sum += c[x + y * res + z * res * res]; + } + } + const real_t rx = dims[0] + (x + 0.5) * bl; + std::cout << rx << "," << sum / (res * res) << "\n"; + } +} + +} // namespace + +int Simulate(int argc, const char** argv) { + Simulation simulation(argc, argv); + + // BDM's Simulation consumes argv and rejects unknown flags, so use + // an env var for the demo-only toggle. Set TRADEWIND_NO_WIND=1 to run + // the pure-diffusion baseline. + const char* flag = std::getenv("TRADEWIND_NO_WIND"); + const bool with_wind = !(flag && std::string(flag) == "1"); + + auto* rm = simulation.GetResourceManager(); + auto* env = simulation.GetEnvironment(); + + // Dummy cells to anchor environment bounds to roughly [0, L]^3. + auto* c0 = new Cell({0, 0, 0}); + c0->SetDiameter(1); + rm->AddAgent(c0); + auto* c1 = new Cell({kL, kL, kL}); + c1->SetDiameter(1); + rm->AddAgent(c1); + env->ForcedUpdate(); + + auto* g = new AdvectionDiffusionGrid(0, "A", kDiff, 0.0, kRes); + g->Initialize(); + g->SetBoundaryConditionType(BoundaryConditionType::kClosedBoundaries); + + SeedRearSource(g); + if (with_wind) SetUniformWind(g, kWind); + InstallBarrier(g); + + PrintXProfile(*g, "t=0"); + std::cout << "# leading_edge_mass,t\n"; + std::cout << "0," << LeadingEdgeMass(*g) << "\n"; + + for (int step = 1; step <= kSteps; step++) { + g->Diffuse(kDt); + if (step % 50 == 0) { + std::cout << step * kDt << "," << LeadingEdgeMass(*g) << "\n"; + } + } + + PrintXProfile(*g, with_wind ? "t_end_wind" : "t_end_nowind"); + + delete g; + return 0; +} + +} // namespace bdm + +int main(int argc, const char** argv) { return bdm::Simulate(argc, argv); } diff --git a/demo/tradewind/src/tradewind.h b/demo/tradewind/src/tradewind.h new file mode 100644 index 000000000..a0c9a99c5 --- /dev/null +++ b/demo/tradewind/src/tradewind.h @@ -0,0 +1,19 @@ +// ----------------------------------------------------------------------------- +// +// Copyright (C) 2026 stanbot8 fork of BioDynaMo. +// Licensed under the Apache License, Version 2.0 (the "License"). +// +// ----------------------------------------------------------------------------- + +#ifndef DEMO_TRADEWIND_H_ +#define DEMO_TRADEWIND_H_ + +#include "biodynamo.h" + +namespace bdm { + +int Simulate(int argc, const char** argv); + +} // namespace bdm + +#endif // DEMO_TRADEWIND_H_