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/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_ diff --git a/src/core/diffusion/advection_diffusion_grid.cc b/src/core/diffusion/advection_diffusion_grid.cc new file mode 100644 index 000000000..420b777bd --- /dev/null +++ b/src/core/diffusion/advection_diffusion_grid.cc @@ -0,0 +1,247 @@ +// ----------------------------------------------------------------------------- +// +// 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(); + + // 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++) { + 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..95d285c1f --- /dev/null +++ b/src/core/diffusion/advection_diffusion_grid.h @@ -0,0 +1,93 @@ +// ----------------------------------------------------------------------------- +// +// 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). +/// +/// 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; + 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(); + + /// 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 + /// 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_ diff --git a/src/core/diffusion/diffusion_grid.h b/src/core/diffusion/diffusion_grid.h index 0b5e74c18..8b22b7543 100644 --- a/src/core/diffusion/diffusion_grid.h +++ b/src/core/diffusion/diffusion_grid.h @@ -348,9 +348,21 @@ 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; 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 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