BioDynaMo  v1.05.120-25dc9790
diffusion_grid.h
Go to the documentation of this file.
1 // -----------------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 CERN & University of Surrey for the benefit of the
4 // BioDynaMo collaboration. All Rights Reserved.
5 //
6 // Licensed under the Apache License, Version 2.0 (the "License");
7 // you may not use this file except in compliance with the License.
8 //
9 // See the LICENSE file distributed with this work for details.
10 // See the NOTICE file distributed with this work for additional information
11 // regarding copyright ownership.
12 //
13 // -----------------------------------------------------------------------------
14 
15 #ifndef CORE_DIFFUSION_DIFFUSION_GRID_H_
16 #define CORE_DIFFUSION_DIFFUSION_GRID_H_
17 
18 #include <array>
19 #include <functional>
20 #include <memory>
21 #include <string>
22 #include <utility>
23 #include <vector>
24 
28 #include "core/util/log.h"
29 #include "core/util/root.h"
30 #include "core/util/spinlock.h"
31 
32 namespace bdm {
33 
36  kDirichlet,
37  kNeumann,
40  kPeriodic
41 };
42 
43 enum class InteractionMode { kAdditive = 0, kExponential = 1, kLogistic = 2 };
44 
50  public:
51  BoundaryCondition() = default;
52  explicit BoundaryCondition(const TRootIOCtor*) {}
53  virtual ~BoundaryCondition() = default;
54 
62  virtual real_t Evaluate(real_t x, real_t y, real_t z, real_t time) const = 0;
63 
65 };
66 
72 
73  public:
76  explicit ConstantBoundaryCondition(real_t value) : value_(value) {}
77 
79  real_t Evaluate(real_t x, real_t y, real_t z, real_t time) const final {
80  return value_;
81  }
82 
83  private:
85  real_t value_ = 0.0;
86 
88 };
89 
90 class DiffusionGrid : public ScalarField {
91  public:
92  DiffusionGrid() = default;
93  explicit DiffusionGrid(const TRootIOCtor*) {}
94  DiffusionGrid(int substance_id, const std::string& substance_name, real_t dc,
95  real_t mu, int resolution = 10);
96  ~DiffusionGrid() override = default;
97  DiffusionGrid(const DiffusionGrid&) = delete; // copy constructor
98  DiffusionGrid& operator=(const DiffusionGrid&) = delete; // copy assignment
99  DiffusionGrid(DiffusionGrid&&) = delete; // move constructor
100  DiffusionGrid& operator=(DiffusionGrid&&) = delete; // move assignment
101 
102  void Initialize() override;
103 
108  void Update() override;
109 
110  void Step(real_t dt) override { Diffuse(dt); }
111  void Diffuse(real_t dt);
112 
113  virtual void DiffuseWithClosedEdge(real_t dt) = 0;
114  virtual void DiffuseWithOpenEdge(real_t dt) = 0;
117  virtual void DiffuseWithDirichlet(real_t dt) = 0;
120  virtual void DiffuseWithNeumann(real_t dt) = 0;
123  virtual void DiffuseWithPeriodic(real_t dt) = 0;
124 
133  void CalculateGradient();
134 
138  void RunInitializers();
139 
157  void ChangeConcentrationBy(const Real3& position, real_t amount,
159  bool scale_with_resolution = false);
162  void ChangeConcentrationBy(size_t idx, real_t amount,
164  bool scale_with_resolution = false);
165 
169  real_t GetValue(const Real3& position) const override;
170  [[deprecated("Use GetValue instead")]] real_t GetConcentration(
171  const Real3& position) const {
172  return GetValue(position);
173  };
177  real_t GetConcentration(const size_t idx) const;
178 
196  Real3 GetGradient(const Real3& position) const override {
197  Real3 gradient;
198  GetGradient(position, &gradient, false);
199  return gradient;
200  };
201 
202  // NOTE: virtual because of test
208  virtual void GetGradient(const Real3& position, Real3* gradient,
209  bool normalize = true) const;
210 
212  std::array<uint32_t, 3> GetBoxCoordinates(const Real3& position) const;
213 
214  // Get the coordinates of the box at the specified index
215  std::array<uint32_t, 3> GetBoxCoordinates(const size_t idx) const;
216 
218  size_t GetBoxIndex(const std::array<uint32_t, 3>& box_coord) const;
219 
221  size_t GetBoxIndex(const Real3& position) const;
222 
227  std::array<size_t, 6> GetNeighboringBoxes(size_t index) const;
228 
233  std::array<size_t, 6> GetNeighboringBoxes(
234  size_t index, const std::array<uint32_t, 3>& box_coord) const;
235 
237  // Check is done in ParametersCheck within the Diffuse(dt) method
238  mu_ = mu;
239  }
240 
242  real_t GetLastTimestep() const { return last_dt_; }
243 
244  // Sets an upper threshold for allowed values in the diffusion grid.
246 
247  // Returns the upper threshold for allowed values in the diffusion grid.
249 
250  // Sets a lower threshold for allowed values in the diffusion grid.
252 
253  // Returns the lower threshold for allowed values in the diffusion grid.
255 
256  const real_t* GetAllConcentrations() const { return c1_.data(); }
257 
258  const real_t* GetAllGradients() const { return gradients_.data()->data(); }
259 
260  std::array<size_t, 3> GetNumBoxesArray() const {
261  std::array<size_t, 3> ret;
262  ret[0] = resolution_;
263  ret[1] = resolution_;
264  ret[2] = resolution_;
265  return ret;
266  }
267 
268  size_t GetNumBoxes() const { return total_num_boxes_; }
269 
270  real_t GetBoxLength() const { return box_length_; }
271 
272  [[deprecated("Use GetContinuumId() instead.")]] int GetSubstanceId() const {
273  return GetContinuumId();
274  }
275 
276  [[deprecated("Use GetContinuumName() instead.")]] const std::string&
278  return GetContinuumName();
279  }
280 
281  real_t GetDecayConstant() const { return mu_; }
282 
283  const int32_t* GetDimensionsPtr() const { return grid_dimensions_.data(); }
284 
285  std::array<int32_t, 6> GetDimensions() const {
286  std::array<int32_t, 6> ret;
287  ret[0] = grid_dimensions_[0];
288  ret[1] = grid_dimensions_[1];
289  ret[2] = grid_dimensions_[0];
290  ret[3] = grid_dimensions_[1];
291  ret[4] = grid_dimensions_[0];
292  ret[5] = grid_dimensions_[1];
293  return ret;
294  }
295 
296  std::array<int32_t, 3> GetGridSize() const {
297  std::array<int32_t, 3> ret;
298  ret[0] = grid_dimensions_[1] - grid_dimensions_[0];
299  ret[1] = grid_dimensions_[1] - grid_dimensions_[0];
300  ret[2] = grid_dimensions_[1] - grid_dimensions_[0];
301  return ret;
302  }
303 
304  const std::array<real_t, 7>& GetDiffusionCoefficients() const { return dc_; }
305 
306  size_t GetResolution() const { return resolution_; }
307 
308  real_t GetBoxVolume() const { return box_volume_; }
309 
310  template <typename F>
311  void AddInitializer(F function) {
312  initializers_.push_back(function);
313  }
314 
317  return (mu_ == 0 && dc_[1] == 0 && dc_[2] == 0 && dc_[3] == 0 &&
318  dc_[4] == 0 && dc_[5] == 0 && dc_[6] == 0);
319  }
320 
324  void SetBoundaryCondition(std::unique_ptr<BoundaryCondition> bc);
325 
329 
332 
335 
337  void PrintInfo(std::ostream& out = std::cout);
338 
342  };
343 
345  bool IsInitialized() const { return initialized_; }
346 
350 
351  private:
352  friend class EulerGrid;
353  friend class EulerDepletionGrid;
354  friend class TestGrid; // class used for testing (e.g. initialization)
355 
356  void ParametersCheck(real_t dt);
357 
370  void CopyOldData(const ParallelResizeVector<real_t>& old_c1,
371  const ParallelResizeVector<Real3>& old_gradients,
372  size_t old_resolution);
373 
392  std::array<real_t, 7> dc_ = {{0}};
394  real_t mu_ = 0;
396  std::array<int32_t, 2> grid_dimensions_ = {{0}};
398  uint64_t num_boxes_axis_ = 0;
400  size_t total_num_boxes_ = 0;
403  size_t resolution_ = 0;
407  bool parity_ = false;
410  std::vector<std::function<real_t(real_t, real_t, real_t)>> initializers_ =
411  {};
412  // Turn to true after gradient initialization
413  bool init_gradient_ = false;
417  std::unique_ptr<BoundaryCondition> boundary_condition_ = nullptr;
419  bool initialized_ = false;
426 
428 };
429 
430 } // namespace bdm
431 
432 #endif // CORE_DIFFUSION_DIFFUSION_GRID_H_
bdm::BoundaryConditionType::kPeriodic
@ kPeriodic
bdm::DiffusionGrid::gradients_
ParallelResizeVector< Real3 > gradients_
The array of gradients (x, y, z)
Definition: diffusion_grid.h:386
bdm::DiffusionGrid::boundary_condition_
std::unique_ptr< BoundaryCondition > boundary_condition_
Object that implements the boundary conditions.
Definition: diffusion_grid.h:417
bdm::DiffusionGrid::Update
void Update() override
Definition: diffusion_grid.cc:156
bdm::DiffusionGrid::GetConcentration
real_t GetConcentration(const Real3 &position) const
Definition: diffusion_grid.h:170
bdm::DiffusionGrid::precompute_gradients_
bool precompute_gradients_
Definition: diffusion_grid.h:425
bdm::DiffusionGrid::GetSubstanceName
const std::string & GetSubstanceName() const
Definition: diffusion_grid.h:277
bdm::DiffusionGrid::SetUpperThreshold
void SetUpperThreshold(real_t t)
Definition: diffusion_grid.h:245
spinlock.h
bdm::DiffusionGrid::PrintInfo
void PrintInfo(std::ostream &out=std::cout)
Print information about the Diffusion Grid.
Definition: diffusion_grid.cc:552
bdm::DiffusionGrid::GetBoundaryCondition
BoundaryCondition * GetBoundaryCondition() const
Returns the boundary condition. Does not transfer ownership.
Definition: diffusion_grid.cc:536
bdm::BoundaryConditionType::kOpenBoundaries
@ kOpenBoundaries
bdm::DiffusionGrid::AddInitializer
void AddInitializer(F function)
Definition: diffusion_grid.h:311
bdm::DiffusionGrid::grid_dimensions_
std::array< int32_t, 2 > grid_dimensions_
The grid dimensions of the diffusion grid (cubic shaped)
Definition: diffusion_grid.h:396
bdm::DiffusionGrid::GetBoxLength
real_t GetBoxLength() const
Definition: diffusion_grid.h:270
bdm::DiffusionGrid::GetValue
real_t GetValue(const Real3 &position) const override
Get the value of the scalar field at specified position.
Definition: diffusion_grid.cc:401
bdm::DiffusionGrid::GetBoxIndex
size_t GetBoxIndex(const std::array< uint32_t, 3 > &box_coord) const
Calculates the box index at specified box coordinates.
Definition: diffusion_grid.cc:497
bdm::DiffusionGrid::parity_
bool parity_
If false, grid dimensions are even; if true, they are odd.
Definition: diffusion_grid.h:407
bdm::DiffusionGrid::CalculateGradient
void CalculateGradient()
Definition: diffusion_grid.cc:317
bdm
Definition: agent.cc:39
bdm::BoundaryCondition::BoundaryCondition
BoundaryCondition()=default
bdm::InteractionMode::kAdditive
@ kAdditive
bdm::DiffusionGrid::GetNumBoxes
size_t GetNumBoxes() const
Definition: diffusion_grid.h:268
bdm::DiffusionGrid::Diffuse
void Diffuse(real_t dt)
Definition: diffusion_grid.cc:125
bdm::DiffusionGrid::print_info_with_initialization_
bool print_info_with_initialization_
Definition: diffusion_grid.h:422
bdm::DiffusionGrid::lower_threshold_
real_t lower_threshold_
The minimum concentration value that a box can have.
Definition: diffusion_grid.h:390
bdm::DiffusionGrid::mu_
real_t mu_
The decay constant.
Definition: diffusion_grid.h:394
bdm::BoundaryConditionType::kClosedBoundaries
@ kClosedBoundaries
bdm::DiffusionGrid::GetLastTimestep
real_t GetLastTimestep() const
Return the last timestep dt that was used to run Diffuse(dt)
Definition: diffusion_grid.h:242
bdm::DiffusionGrid::locks_
ParallelResizeVector< Spinlock > locks_
Definition: diffusion_grid.h:380
bdm::real_t
double real_t
Definition: real_t.h:21
bdm::DiffusionGrid::GetBoxVolume
real_t GetBoxVolume() const
Definition: diffusion_grid.h:308
bdm::DiffusionGrid::BDM_CLASS_DEF_OVERRIDE
BDM_CLASS_DEF_OVERRIDE(DiffusionGrid, 1)
bdm::DiffusionGrid::DiffuseWithOpenEdge
virtual void DiffuseWithOpenEdge(real_t dt)=0
bdm::DiffusionGrid::GetDecayConstant
real_t GetDecayConstant() const
Definition: diffusion_grid.h:281
bdm::DiffusionGrid::c1_
ParallelResizeVector< real_t > c1_
The array of concentration values.
Definition: diffusion_grid.h:382
bdm::DiffusionGrid::box_volume_
real_t box_volume_
the volume of each box
Definition: diffusion_grid.h:377
bdm::InteractionMode::kLogistic
@ kLogistic
bdm::DiffusionGrid::init_gradient_
bool init_gradient_
Definition: diffusion_grid.h:413
bdm::DiffusionGrid::GetBoxCoordinates
std::array< uint32_t, 3 > GetBoxCoordinates(const Real3 &position) const
Get the coordinates of the box at the specified position.
Definition: diffusion_grid.cc:463
bdm::DiffusionGrid::SetBoundaryConditionType
void SetBoundaryConditionType(BoundaryConditionType bc_type)
Sets boundary condition type.
Definition: diffusion_grid.cc:540
bdm::ParallelResizeVector< real_t >
bdm::DiffusionGrid::ParametersCheck
void ParametersCheck(real_t dt)
Definition: diffusion_grid.cc:593
bdm::DiffusionGrid::total_num_boxes_
size_t total_num_boxes_
The total number of boxes in the diffusion grid.
Definition: diffusion_grid.h:400
bdm::DiffusionGrid::GetUpperThreshold
real_t GetUpperThreshold() const
Definition: diffusion_grid.h:248
bdm::DiffusionGrid::CopyOldData
void CopyOldData(const ParallelResizeVector< real_t > &old_c1, const ParallelResizeVector< Real3 > &old_gradients, size_t old_resolution)
Definition: diffusion_grid.cc:211
bdm::DiffusionGrid
Definition: diffusion_grid.h:90
bdm::DiffusionGrid::GetGridSize
std::array< int32_t, 3 > GetGridSize() const
Definition: diffusion_grid.h:296
bdm::BoundaryCondition::BoundaryCondition
BoundaryCondition(const TRootIOCtor *)
Definition: diffusion_grid.h:52
bdm::BoundaryConditionType
BoundaryConditionType
Available boundary conditions.
Definition: diffusion_grid.h:35
bdm::InteractionMode
InteractionMode
Definition: diffusion_grid.h:43
bdm::DiffusionGrid::GetSubstanceId
int GetSubstanceId() const
Definition: diffusion_grid.h:272
bdm::Continuum::GetContinuumName
const std::string & GetContinuumName() const
Returns the name of the continuum.
Definition: continuum_interface.h:93
bdm::DiffusionGrid::DiffuseWithPeriodic
virtual void DiffuseWithPeriodic(real_t dt)=0
bdm::BoundaryCondition::Evaluate
virtual real_t Evaluate(real_t x, real_t y, real_t z, real_t time) const =0
Boundary condition for Neumann and Dirichlet boundary conditions.
bdm::DiffusionGrid::DiffusionGrid
DiffusionGrid(const TRootIOCtor *)
Definition: diffusion_grid.h:93
bdm::DiffusionGrid::DiffuseWithClosedEdge
virtual void DiffuseWithClosedEdge(real_t dt)=0
bdm::ConstantBoundaryCondition::Evaluate
real_t Evaluate(real_t x, real_t y, real_t z, real_t time) const final
see BoundaryCondition::Evaluate()
Definition: diffusion_grid.h:79
parallel_resize_vector.h
bdm::DiffusionGrid::GetBoundaryConditionType
BoundaryConditionType GetBoundaryConditionType() const
Returns the BoundaryConditionType, see BoundaryConditionType
Definition: diffusion_grid.cc:548
bdm::ConstantBoundaryCondition::ConstantBoundaryCondition
ConstantBoundaryCondition(real_t value)
Constructor for a constant boundary condition.
Definition: diffusion_grid.h:76
bdm::DiffusionGrid::~DiffusionGrid
~DiffusionGrid() override=default
bdm::InteractionMode::kExponential
@ kExponential
bdm::DiffusionGrid::RunInitializers
void RunInitializers()
Definition: diffusion_grid.cc:253
continuum_interface.h
bdm::BoundaryCondition
This class implements the boundary conditions. It is integrated into the diffusion grid as a smart po...
Definition: diffusion_grid.h:49
bdm::DiffusionGrid::initializers_
std::vector< std::function< real_t(real_t, real_t, real_t)> > initializers_
Definition: diffusion_grid.h:410
bdm::DiffusionGrid::initialized_
bool initialized_
Flag to indicate if the grid is initialized.
Definition: diffusion_grid.h:419
bdm::Continuum::GetContinuumId
int GetContinuumId() const
Returns the ID of the continuum.
Definition: continuum_interface.h:87
root.h
bdm::EulerDepletionGrid
Continuum model for the 3D diffusion equation with exponential decay and substance depletion .
Definition: euler_depletion_grid.h:34
bdm::DiffusionGrid::Initialize
void Initialize() override
Definition: diffusion_grid.cc:77
bdm::DiffusionGrid::GetLowerThreshold
real_t GetLowerThreshold() const
Definition: diffusion_grid.h:254
math_array.h
bdm::DiffusionGrid::Step
void Step(real_t dt) override
Definition: diffusion_grid.h:110
bdm::DiffusionGrid::GetNeighboringBoxes
std::array< size_t, 6 > GetNeighboringBoxes(size_t index) const
Definition: diffusion_grid.cc:510
bdm::ParallelResizeVector::data
T * data() noexcept
Definition: parallel_resize_vector.h:74
bdm::DiffusionGrid::box_length_
real_t box_length_
The side length of each box.
Definition: diffusion_grid.h:375
bdm::DiffusionGrid::PrintInfoWithInitialization
void PrintInfoWithInitialization()
Print the information after initialization.
Definition: diffusion_grid.h:340
bdm::DiffusionGrid::GetNumBoxesArray
std::array< size_t, 3 > GetNumBoxesArray() const
Definition: diffusion_grid.h:260
bdm::DiffusionGrid::GetDiffusionCoefficients
const std::array< real_t, 7 > & GetDiffusionCoefficients() const
Definition: diffusion_grid.h:304
bdm::ScalarField
Interface for scalar fields. See Continuum for more information.
Definition: continuum_interface.h:132
bdm::BoundaryCondition::~BoundaryCondition
virtual ~BoundaryCondition()=default
log.h
bdm::DiffusionGrid::TestGrid
friend class TestGrid
Definition: diffusion_grid.h:354
bdm::DiffusionGrid::DiffusionGrid
DiffusionGrid()=default
bdm::DiffusionGrid::last_dt_
real_t last_dt_
The last timestep dt used for the diffusion grid update Diffuse(dt)
Definition: diffusion_grid.h:405
bdm::ConstantBoundaryCondition::BDM_CLASS_DEF_OVERRIDE
BDM_CLASS_DEF_OVERRIDE(ConstantBoundaryCondition, 1)
bdm::DiffusionGrid::GetDimensions
std::array< int32_t, 6 > GetDimensions() const
Definition: diffusion_grid.h:285
bdm::DiffusionGrid::bc_type_
BoundaryConditionType bc_type_
Type of boundary conditions.
Definition: diffusion_grid.h:415
bdm::BoundaryConditionType::kNeumann
@ kNeumann
bdm::DiffusionGrid::resolution_
size_t resolution_
Definition: diffusion_grid.h:403
bdm::DiffusionGrid::GetAllConcentrations
const real_t * GetAllConcentrations() const
Definition: diffusion_grid.h:256
bdm::DiffusionGrid::TurnOffGradientCalculation
void TurnOffGradientCalculation()
Definition: diffusion_grid.h:349
bdm::DiffusionGrid::SetDecayConstant
void SetDecayConstant(real_t mu)
Definition: diffusion_grid.h:236
bdm::DiffusionGrid::SetLowerThreshold
void SetLowerThreshold(real_t t)
Definition: diffusion_grid.h:251
bdm::DiffusionGrid::DiffuseWithDirichlet
virtual void DiffuseWithDirichlet(real_t dt)=0
bdm::DiffusionGrid::DiffuseWithNeumann
virtual void DiffuseWithNeumann(real_t dt)=0
bdm::DiffusionGrid::GetDimensionsPtr
const int32_t * GetDimensionsPtr() const
Definition: diffusion_grid.h:283
bdm::MathArray
Definition: math_array.h:36
bdm::BoundaryCondition::BDM_CLASS_DEF
BDM_CLASS_DEF(BoundaryCondition, 1)
bdm::DiffusionGrid::upper_threshold_
real_t upper_threshold_
The maximum concentration value that a box can have.
Definition: diffusion_grid.h:388
bdm::DiffusionGrid::IsInitialized
bool IsInitialized() const
Returns if the grid has been initialized.
Definition: diffusion_grid.h:345
bdm::DiffusionGrid::ChangeConcentrationBy
void ChangeConcentrationBy(const Real3 &position, real_t amount, InteractionMode mode=InteractionMode::kAdditive, bool scale_with_resolution=false)
Definition: diffusion_grid.cc:355
bdm::DiffusionGrid::SetBoundaryCondition
void SetBoundaryCondition(std::unique_ptr< BoundaryCondition > bc)
Set the boundary condition, takes ownership of the object.
Definition: diffusion_grid.cc:531
bdm::BoundaryConditionType::kDirichlet
@ kDirichlet
bdm::DiffusionGrid::operator=
DiffusionGrid & operator=(const DiffusionGrid &)=delete
bdm::DiffusionGrid::GetGradient
Real3 GetGradient(const Real3 &position) const override
Definition: diffusion_grid.h:196
bdm::ConstantBoundaryCondition
This class implements constant boundary conditions (Dirichlet and Neumann). The value of the boundary...
Definition: diffusion_grid.h:70
bdm::DiffusionGrid::dc_
std::array< real_t, 7 > dc_
The diffusion coefficients [cc, cw, ce, cs, cn, cb, ct].
Definition: diffusion_grid.h:392
bdm::DiffusionGrid::GetAllGradients
const real_t * GetAllGradients() const
Definition: diffusion_grid.h:258
bdm::DiffusionGrid::c2_
ParallelResizeVector< real_t > c2_
An extra concentration data buffer for faster value updating.
Definition: diffusion_grid.h:384
bdm::DiffusionGrid::GetResolution
size_t GetResolution() const
Definition: diffusion_grid.h:306
bdm::DiffusionGrid::num_boxes_axis_
uint64_t num_boxes_axis_
The number of boxes at each axis [x, y, z] (same along each axis)
Definition: diffusion_grid.h:398
bdm::ConstantBoundaryCondition::value_
real_t value_
Constant value of the boundary condition for all positions and times.
Definition: diffusion_grid.h:85
bdm::DiffusionGrid::IsFixedSubstance
bool IsFixedSubstance()
Return if a substance is stationary, e.g. (mu == 0 && dc == 0)
Definition: diffusion_grid.h:316
bdm::EulerGrid
Continuum model for the 3D heat equation with exponential decay .
Definition: euler_grid.h:43