diffusion_grid.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// Copyright (C) 2021 CERN & University of Surrey for the benefit of the
// BioDynaMo collaboration. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
//
// See the LICENSE file distributed with this work for details.
// See the NOTICE file distributed with this work for additional information
// regarding copyright ownership.
//
// -----------------------------------------------------------------------------
 
#ifndef CORE_DIFFUSION_DIFFUSION_GRID_H_
#define CORE_DIFFUSION_DIFFUSION_GRID_H_
 
#include <array>
#include <functional>
#include <memory>
#include <string>
#include <utility>
#include <vector>
 
#include "core/container/math_array.h"
#include "core/container/parallel_resize_vector.h"
#include "core/diffusion/continuum_interface.h"
#include "core/util/log.h"
#include "core/util/root.h"
#include "core/util/spinlock.h"
 
namespace bdm {
 
enum class BoundaryConditionType {
  kDirichlet,
  kNeumann,
  kOpenBoundaries,
  kClosedBoundaries,
  kPeriodic
};
 
enum class InteractionMode { kAdditive = 0, kExponential = 1, kLogistic = 2 };
 
class BoundaryCondition {
 public:
  BoundaryCondition() = default;
  explicit BoundaryCondition(const TRootIOCtor*) {}
  virtual ~BoundaryCondition() = default;
 
  virtual real_t Evaluate(real_t x, real_t y, real_t z, real_t time) const = 0;
 
  BDM_CLASS_DEF(BoundaryCondition, 1);
};
 
class ConstantBoundaryCondition : public BoundaryCondition {
  using BoundaryCondition::BoundaryCondition;
 
 public:
  explicit ConstantBoundaryCondition(real_t value) : value_(value) {}
 
  real_t Evaluate(real_t x, real_t y, real_t z, real_t time) const final {
    return value_;
  }
 
 private:
  real_t value_ = 0.0;
 
  BDM_CLASS_DEF_OVERRIDE(ConstantBoundaryCondition, 1);
};
 
class DiffusionGrid : public ScalarField {
 public:
  DiffusionGrid() = default;
  explicit DiffusionGrid(const TRootIOCtor*) {}
  DiffusionGrid(int substance_id, const std::string& substance_name, real_t dc,
                real_t mu, int resolution = 10);
  ~DiffusionGrid() override = default;
  DiffusionGrid(const DiffusionGrid&) = delete;             // copy constructor
  DiffusionGrid& operator=(const DiffusionGrid&) = delete;  // copy assignment
  DiffusionGrid(DiffusionGrid&&) = delete;                  // move constructor
  DiffusionGrid& operator=(DiffusionGrid&&) = delete;       // move assignment
 
  void Initialize() override;
 
  void Update() override;
 
  void Step(real_t dt) override { Diffuse(dt); }
  void Diffuse(real_t dt);
 
  virtual void DiffuseWithClosedEdge(real_t dt) = 0;
  virtual void DiffuseWithOpenEdge(real_t dt) = 0;
  virtual void DiffuseWithDirichlet(real_t dt) = 0;
  virtual void DiffuseWithNeumann(real_t dt) = 0;
  virtual void DiffuseWithPeriodic(real_t dt) = 0;
 
  void CalculateGradient();
 
  void RunInitializers();
 
  void ChangeConcentrationBy(const Real3& position, real_t amount,
                             InteractionMode mode = InteractionMode::kAdditive,
                             bool scale_with_resolution = false);
  void ChangeConcentrationBy(size_t idx, real_t amount,
                             InteractionMode mode = InteractionMode::kAdditive,
                             bool scale_with_resolution = false);
 
  real_t GetValue(const Real3& position) const override;
  [[deprecated("Use GetValue instead")]] real_t GetConcentration(
      const Real3& position) const {
    return GetValue(position);
  };
  real_t GetConcentration(const size_t idx) const;
 
  Real3 GetGradient(const Real3& position) const override {
    Real3 gradient;
    GetGradient(position, &gradient, false);
    return gradient;
  };
 
  // NOTE: virtual because of test
  virtual void GetGradient(const Real3& position, Real3* gradient,
                           bool normalize = true) const;
 
  std::array<uint32_t, 3> GetBoxCoordinates(const Real3& position) const;
 
  // Get the coordinates of the box at the specified index
  std::array<uint32_t, 3> GetBoxCoordinates(const size_t idx) const;
 
  size_t GetBoxIndex(const std::array<uint32_t, 3>& box_coord) const;
 
  size_t GetBoxIndex(const Real3& position) const;
 
  std::array<size_t, 6> GetNeighboringBoxes(size_t index) const;
 
  std::array<size_t, 6> GetNeighboringBoxes(
      size_t index, const std::array<uint32_t, 3>& box_coord) const;
 
  void SetDecayConstant(real_t mu) {
    // Check is done in ParametersCheck within the Diffuse(dt) method
    mu_ = mu;
  }
 
  real_t GetLastTimestep() const { return last_dt_; }
 
  // Sets an upper threshold for allowed values in the diffusion grid.
  void SetUpperThreshold(real_t t) { upper_threshold_ = t; }
 
  // Returns the upper threshold for allowed values in the diffusion grid.
  real_t GetUpperThreshold() const { return upper_threshold_; }
 
  // Sets a lower threshold for allowed values in the diffusion grid.
  void SetLowerThreshold(real_t t) { lower_threshold_ = t; }
 
  // Returns the lower threshold for allowed values in the diffusion grid.
  real_t GetLowerThreshold() const { return lower_threshold_; }
 
  const real_t* GetAllConcentrations() const { return c1_.data(); }
 
  const real_t* GetAllGradients() const { return gradients_.data()->data(); }
 
  std::array<size_t, 3> GetNumBoxesArray() const {
    std::array<size_t, 3> ret;
    ret[0] = resolution_;
    ret[1] = resolution_;
    ret[2] = resolution_;
    return ret;
  }
 
  size_t GetNumBoxes() const { return total_num_boxes_; }
 
  real_t GetBoxLength() const { return box_length_; }
 
  [[deprecated("Use GetContinuumId() instead.")]] int GetSubstanceId() const {
    return GetContinuumId();
  }
 
  [[deprecated("Use GetContinuumName() instead.")]] const std::string&
  GetSubstanceName() const {
    return GetContinuumName();
  }
 
  real_t GetDecayConstant() const { return mu_; }
 
  const int32_t* GetDimensionsPtr() const { return grid_dimensions_.data(); }
 
  std::array<int32_t, 6> GetDimensions() const {
    std::array<int32_t, 6> ret;
    ret[0] = grid_dimensions_[0];
    ret[1] = grid_dimensions_[1];
    ret[2] = grid_dimensions_[0];
    ret[3] = grid_dimensions_[1];
    ret[4] = grid_dimensions_[0];
    ret[5] = grid_dimensions_[1];
    return ret;
  }
 
  std::array<int32_t, 3> GetGridSize() const {
    std::array<int32_t, 3> ret;
    ret[0] = grid_dimensions_[1] - grid_dimensions_[0];
    ret[1] = grid_dimensions_[1] - grid_dimensions_[0];
    ret[2] = grid_dimensions_[1] - grid_dimensions_[0];
    return ret;
  }
 
  const std::array<real_t, 7>& GetDiffusionCoefficients() const { return dc_; }
 
  size_t GetResolution() const { return resolution_; }
 
  real_t GetBoxVolume() const { return box_volume_; }
 
  template <typename F>
  void AddInitializer(F function) {
    initializers_.push_back(function);
  }
 
  bool IsFixedSubstance() {
    return (mu_ == 0 && dc_[1] == 0 && dc_[2] == 0 && dc_[3] == 0 &&
            dc_[4] == 0 && dc_[5] == 0 && dc_[6] == 0);
  }
 
  void SetBoundaryCondition(std::unique_ptr<BoundaryCondition> bc);
 
  BoundaryCondition* GetBoundaryCondition() const;
 
  void SetBoundaryConditionType(BoundaryConditionType bc_type);
 
  BoundaryConditionType GetBoundaryConditionType() const;
 
  void PrintInfo(std::ostream& out = std::cout);
 
  void PrintInfoWithInitialization() {
    print_info_with_initialization_ = true;
  };
 
  bool IsInitialized() const { return initialized_; }
 
  void TurnOffGradientCalculation() { precompute_gradients_ = false; }
 
 private:
  friend class EulerGrid;
  friend class EulerDepletionGrid;
  friend class TestGrid;  // class used for testing (e.g. initialization)
 
  void ParametersCheck(real_t dt);
 
  void CopyOldData(const ParallelResizeVector<real_t>& old_c1,
                   const ParallelResizeVector<Real3>& old_gradients,
                   size_t old_resolution);
 
  real_t box_length_ = 0;
  real_t box_volume_ = 0;
  mutable ParallelResizeVector<Spinlock> locks_ = {};  
  ParallelResizeVector<real_t> c1_ = {};
  ParallelResizeVector<real_t> c2_ = {};
  ParallelResizeVector<Real3> gradients_ = {};
  real_t upper_threshold_ = 1e15;
  real_t lower_threshold_ = 0.0;
  std::array<real_t, 7> dc_ = {{0}};
  real_t mu_ = 0;
  std::array<int32_t, 2> grid_dimensions_ = {{0}};
  uint64_t num_boxes_axis_ = 0;
  size_t total_num_boxes_ = 0;
  size_t resolution_ = 0;
  real_t last_dt_ = 0.0;
  bool parity_ = false;
  std::vector<std::function<real_t(real_t, real_t, real_t)>> initializers_ =
      {};  
  // Turn to true after gradient initialization
  bool init_gradient_ = false;
  BoundaryConditionType bc_type_ = BoundaryConditionType::kDirichlet;
  std::unique_ptr<BoundaryCondition> boundary_condition_ = nullptr;
  bool initialized_ = false;
  bool print_info_with_initialization_ = false;
  bool precompute_gradients_ = true;
 
  BDM_CLASS_DEF_OVERRIDE(DiffusionGrid, 1);
};
 
}  // namespace bdm
 
#endif  // CORE_DIFFUSION_DIFFUSION_GRID_H_