Skip to content

Skip FTCS solver for zero-diffusion grids #459

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

Open
stanbot8 wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Skip FTCS solver for zero-diffusion grids #459

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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: 12 additions & 4 deletions src/core/diffusion/diffusion_grid.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ class DiffusionGrid : public ScalarField {
[[deprecated("Use GetValue instead")]] real_t GetConcentration(
const Real3& position) const {
return GetValue(position);
};
}
/// @brief Get the concentration at specified index
/// @param idx Flat index of the grid
/// @return c1_[idx]
Expand Down Expand Up @@ -337,9 +337,7 @@ class DiffusionGrid : public ScalarField {
void PrintInfo(std::ostream& out = std::cout);

/// Print the information after initialization
void PrintInfoWithInitialization() {
print_info_with_initialization_ = true;
};
void PrintInfoWithInitialization() { print_info_with_initialization_ = true; }

/// Returns if the grid has been initialized
bool IsInitialized() const { return initialized_; }
Expand Down Expand Up @@ -371,6 +369,16 @@ class DiffusionGrid : public ScalarField {
const ParallelResizeVector<Real3>& old_gradients,
size_t old_resolution);

/// Apply decay without diffusion (used when diffusion coefficient is 0).
void ApplyDecayOnly(real_t dt) {
const real_t decay = 1 - mu_ * dt;
#pragma omp parallel for
for (size_t i = 0; i < total_num_boxes_; i++) {
c2_[i] = c1_[i] * decay;
}
c1_.swap(c2_);
}

/// The side length of each box
real_t box_length_ = 0;
/// the volume of each box
Expand Down
35 changes: 30 additions & 5 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,17 @@
namespace bdm {

void EulerGrid::DiffuseWithClosedEdge(real_t dt) {
const real_t d = 1 - dc_[0];
if (d == 0) {
ApplyDecayOnly(dt);
return;
}

const auto nx = resolution_;
const auto ny = resolution_;
const auto nz = resolution_;

const real_t ibl2 = 1 / (box_length_ * box_length_);
const real_t d = 1 - dc_[0];

constexpr size_t YBF = 16;
#pragma omp parallel for collapse(2)
Expand Down Expand Up @@ -67,12 +72,17 @@ void EulerGrid::DiffuseWithClosedEdge(real_t dt) {
}

void EulerGrid::DiffuseWithOpenEdge(real_t dt) {
const real_t d = 1 - dc_[0];
if (d == 0) {
ApplyDecayOnly(dt);
return;
}

const auto nx = resolution_;
const auto ny = resolution_;
const auto nz = resolution_;

const real_t ibl2 = 1 / (box_length_ * box_length_);
const real_t d = 1 - dc_[0];
std::array<int, 4> l;

constexpr size_t YBF = 16;
Expand Down Expand Up @@ -155,12 +165,17 @@ void EulerGrid::DiffuseWithOpenEdge(real_t dt) {
}

void EulerGrid::DiffuseWithDirichlet(real_t dt) {
const real_t d = 1 - dc_[0];
if (d == 0) {
ApplyDecayOnly(dt);
return;
}

const auto nx = resolution_;
const auto ny = resolution_;
const auto nz = resolution_;

const real_t ibl2 = 1 / (box_length_ * box_length_);
const real_t d = 1 - dc_[0];

const auto sim_time = GetSimulatedTime();

Expand Down Expand Up @@ -211,13 +226,18 @@ void EulerGrid::DiffuseWithDirichlet(real_t dt) {
}

void EulerGrid::DiffuseWithNeumann(real_t dt) {
const real_t d = 1 - dc_[0];
if (d == 0) {
ApplyDecayOnly(dt);
return;
}

const size_t nx = resolution_;
const size_t ny = resolution_;
const size_t nz = resolution_;
const size_t num_boxes = nx * ny * nz;

const real_t ibl2 = 1 / (box_length_ * box_length_);
const real_t d = 1 - dc_[0];

const auto sim_time = GetSimulatedTime();

Expand Down Expand Up @@ -302,12 +322,17 @@ void EulerGrid::DiffuseWithNeumann(real_t dt) {
}

void EulerGrid::DiffuseWithPeriodic(real_t dt) {
const real_t d = 1 - dc_[0];
if (d == 0) {
ApplyDecayOnly(dt);
return;
}

const size_t nx = resolution_;
const size_t ny = resolution_;
const size_t nz = resolution_;

const real_t dx = box_length_;
const real_t d = 1 - dc_[0];

constexpr size_t YBF = 16;
#pragma omp parallel for collapse(2)
Expand Down
Loading