monod.hpp

#ifndef __MODEL_MONOD_HPP__
#define __MODEL_MONOD_HPP__
 
#include "mc/macros.hpp"
#include "models/utils.hpp"
#include <Kokkos_MathematicalConstants.hpp>
#include <mc/alias.hpp>
#include <mc/prng/prng_extension.hpp>
#include <mc/traits.hpp>
#include <optional>
 
namespace Models
{

  struct Monod
  {
    enum class particle_var : int
    {
      l = 0, 
      l_max, 
      mu_p,  
      mue,   
      _init_only_cell_lenghtening, 
      phi_s_c,                     
      __COUNT__ 
    };
    // Monod model consumes only glucose
    MODEL_CONSTANT std::size_t n_c = 1;
    // clang-format off
    static constexpr std::size_t n_var = INDEX_FROM_ENUM(particle_var::__COUNT__);
    static constexpr std::string_view name = "monod";
    using uniform_weight = std::true_type; // Using type alias
    using Self = Monod;
    using FloatType = float;
 
    using SelfParticle = MC::ParticlesModel<Self::n_var, Self::FloatType>;
    using SelfContribs = MC::ParticlesContribs<Self::n_c, Self::FloatType>;
    using Config = std::nullopt_t;
 

    MODEL_CONSTANT FloatType y_s_x = 2;                
    MODEL_CONSTANT FloatType mu_max = 0.77 / 3600.;    
    MODEL_CONSTANT FloatType tau_meta = 1. / mu_max;   
    MODEL_CONSTANT FloatType l_max_m = 2e-6;           
    MODEL_CONSTANT FloatType l_min_m = l_max_m / 2.;   
    MODEL_CONSTANT FloatType k_s = 1e-3;               
    MODEL_CONSTANT FloatType d_m = 0.6e-6;             
    MODEL_CONSTANT FloatType lin_density =
        c_linear_density(static_cast<FloatType>(1000), d_m); 
    // clang-format on
 
    MODEL_CONSTANT auto initial_length_dist
        = MC::Distributions::TruncatedNormal<FloatType>(
            l_max_m * 0.75, l_max_m * 0.75 / 4, l_min_m, l_max_m);
 
    KOKKOS_INLINE_FUNCTION static void
    init([[maybe_unused]] const MC::pool_type& random_pool,
         [[maybe_unused]] std::size_t idx,
         [[maybe_unused]] const SelfParticle& arr)
    {
      // Helper cause, initial_length_dist may not be captured in cuda context
      // (FIXME)
      static constexpr auto local_dist = initial_length_dist;
 
      const auto [l0, mu] = MC::sample_random_variables(
          random_pool,
          [](auto& gen)
          {
            const auto l0 = local_dist.draw(
                gen); // Get initial fro given distribution (normal)
            const auto mu = mu_max; // Set maximal for each cell mu_p
            return std::make_tuple(l0, mu);
          });
 
      GET_PROPERTY(particle_var::l) = l0;
      GET_PROPERTY(particle_var::l_max)
          = l_max_m; // Set the same length for everyone
      GET_PROPERTY(particle_var::mu_p) = mu_max;
 
      // born at l_max/2, divides at l_max
      // Need to be changed if division criteria is changed
      constexpr auto dl = l_max_m / 2.;
 
      // Monod factor to convert growth rate into elongation, it assumes
      // doubling in mass during generation time G=ln(2)/mu factor= DL/ln(2)
      // comes from mu=ln(2)/G
      GET_PROPERTY(particle_var::_init_only_cell_lenghtening)
          = dl / Kokkos::numbers::ln2;
    }
 
    KOKKOS_INLINE_FUNCTION static double
    mass([[maybe_unused]] std::size_t idx,
         [[maybe_unused]] const SelfParticle& arr)
    {
      return GET_PROPERTY(Self::particle_var::l) * lin_density; // Rod-shaped
    }
 
    KOKKOS_INLINE_FUNCTION static MC::Status
    update([[maybe_unused]] const MC::pool_type& random_pool,
           FloatType d_t,
           std::size_t idx,
           const SelfParticle& arr,
           const SelfContribs& arr_contribs,
           const std::size_t position_index,
           const MC::LocalConcentration& c)
    {
      constexpr FloatType inv_tau = 1.0 / tau_meta;
 
      const auto s = GET_CLAMPED_CONCENTRATION_CAST(FloatType, 0);
 
      auto& l = GET_PROPERTY(Self::particle_var::l);
      auto& mu_p = GET_PROPERTY(Self::particle_var::mu_p);
 
      // Instantaneous mu from Monod
      const FloatType mu = mu_max * s / (k_s + s);
 
      // Efffective growth rate
      const FloatType mu_eff = Kokkos::min(mu_p, mu);
      // Store only for being exported in 1/h
      GET_PROPERTY(Self::particle_var::mue) = mu_eff * 3600.;
 
      // Contributions
      const auto phi_s
          = -mu_eff * y_s_x * static_cast<FloatType>(mass(idx, arr));
 
      GET_PROPERTY(Self::particle_var::phi_s_c) = phi_s;
      GET_CONTRIBS(0) = phi_s;
 
      // ODE
 
      // Lengthening
      l += d_t
           * (mu_eff
              * GET_PROPERTY(Self::particle_var::_init_only_cell_lenghtening));
 
      // Growth rate
      mu_p += d_t * inv_tau * (mu - mu_p);
 
      return check_div(l, GET_PROPERTY(Self::particle_var::l_max));
    }
 
    KOKKOS_INLINE_FUNCTION static void
    division([[maybe_unused]] const MC::pool_type& random_pool,
             [[maybe_unused]] std::size_t idx,
             [[maybe_unused]] std::size_t idx2,
             [[maybe_unused]] const SelfParticle& arr,
             [[maybe_unused]] const SelfParticle& buffer_arr)
    {
 
      //_init_only_cell_lenghtening doesnt need to be redistributed bcause
      // division criteria is the same for each cell then dl remains the same
 
      // Even if division is deterministic, keep generate way to divide

 
      const FloatType current_l = GET_PROPERTY(particle_var::l);
      const FloatType new_current_length = current_l / 2.F;
      GET_PROPERTY(particle_var::l) = new_current_length;
 
      const FloatType mu_p_value = GET_PROPERTY(particle_var::mu_p);
      const FloatType cell_lenghtening_value
          = GET_PROPERTY(particle_var::_init_only_cell_lenghtening);
 
      GET_PROPERTY_FROM(idx2, buffer_arr, particle_var::l) = new_current_length;
      GET_PROPERTY_FROM(idx2, buffer_arr, particle_var::l_max) = l_max_m;
      GET_PROPERTY_FROM(idx2, buffer_arr, particle_var::mu_p) = mu_p_value;
      GET_PROPERTY_FROM(
          idx2, buffer_arr, particle_var::_init_only_cell_lenghtening)
          = cell_lenghtening_value;
    }
 
    static std::vector<std::string_view>
    names()
    {
      return { "length", "mu", "mu_eff" };
    }
 
    static std::vector<std::size_t>
    get_number()
    {
      return { INDEX_FROM_ENUM(particle_var::l),
               INDEX_FROM_ENUM(particle_var::mu_p),
               INDEX_FROM_ENUM(particle_var::mue) };
    }
 
    static std::vector<std::string_view>
    species()
    {
      return { "S" };
    }
  };
 
  CHECK_MODEL(Monod)
} // namespace Models
#endif