Skip to content

Expose DiffusionGrid buffers via protected accessors; migrate Euler grids off friend pattern #7

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
stanbot8 wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Expose DiffusionGrid buffers via protected accessors; migrate Euler grids off friend pattern #7

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f00c815
Expose DiffusionGrid buffer access to subclasses via protected accessors
stanbot8 Apr 13, 2026
417067b
Convert EulerGrid and EulerDepletionGrid to use protected accessors
stanbot8 Apr 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 14 additions & 2 deletions src/core/diffusion/diffusion_grid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
13 changes: 8 additions & 5 deletions src/core/diffusion/euler_depletion_grid.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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();

Expand All @@ -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) {
Expand All @@ -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) {
Expand Down
186 changes: 108 additions & 78 deletions src/core/diffusion/euler_grid.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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<int, 4> l;

constexpr size_t YBF = 16;
Expand Down Expand Up @@ -122,45 +134,51 @@ 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;
++n;
++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();

Expand All @@ -185,39 +203,46 @@ 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;
s = c + nx;
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();

Expand Down Expand Up @@ -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;
Expand All @@ -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)
Expand Down Expand Up @@ -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
Loading