two_meta_nb.hpp

// #ifndef __TWOMETA_MODEL_NB_HPP__
// #define __TWOMETA_MODEL_NB_HPP__
 
// #include <common/traits.hpp>
// #include <mc/macros.hpp>
// #include <mc/prng/prng_extension.hpp>
// #include <mc/traits.hpp>
// #include <models/uptake_dyn.hpp>
// #include <models/utils.hpp>
// #include <string_view>
 
// namespace Models
// {
 
//   struct TwoMetaNb
//   {
//     using uniform_weight = std::true_type; // Using type alias
//     using Self = TwoMetaNb;
//     using FloatType = float;
//     using Config = std::nullopt_t;
 
//     enum class particle_var : int
//     {
//       age = INDEX_FROM_ENUM(Uptakeparticle_var::COUNT),
//       length,
//       nu1,
//       nu2,
//       l_cp,
//       nu_eff_1, // This is not itself a model property but stored to be
//       exported nu_eff_2, // This is not itself a model property but stored to
//       be exported
//       // TODO FIND BETTER WAY TO STORE/GET CONTRIBUTIONS
//       phi_pts,
//       contrib_phi_s,
//       contrib_phi_o2,
//       contrib_phi_ac,
//       COUNT
//     };
 
//     static constexpr std::size_t n_var =
//     INDEX_FROM_ENUM(particle_var::COUNT); static constexpr std::string_view
//     name = "two_mode_nb"; using SelfParticle =
//     MC::ParticlesModel<Self::n_var, Self::FloatType>;
 
//     // Constants BEGIN
//     MODEL_CONSTANT FloatType MolarMassG
//         = Models::MolarMass::GramPerMole::glucose<FloatType>;
//     MODEL_CONSTANT FloatType MolarMassO2
//         = Models::MolarMass::GramPerMole::dioxygen<FloatType>; // g/mol
 
//     MODEL_CONSTANT FloatType l_max_m = 5e-6;   // m
//     MODEL_CONSTANT FloatType l_c_m = 3e-6;     // m
//     MODEL_CONSTANT FloatType d_m = 0.3e-6;     // m
//     MODEL_CONSTANT FloatType l_min_m = 0.9e-6; // m
 
//     MODEL_CONSTANT FloatType dl_max_ms
//         = 10 * 2e-10; // m/s  https://doi.org/10.7554/eLife.67495;
 
//     MODEL_CONSTANT FloatType lin_density
//         = c_linear_density(static_cast<FloatType>(1000), d_m);
 
//     MODEL_CONSTANT FloatType y_sx_1
//         = 1. / 2.217737e+00;                       // Mode 1 S to X yield
//         (mass)
//     MODEL_CONSTANT FloatType y_sx_2 = y_sx_1 / 3.; // Mode 2 S to X yield
//     (mass) MODEL_CONSTANT FloatType y_sa = 0.8;           // S to A yield
//     (mass) MODEL_CONSTANT FloatType y_os_molar = 3; // 3 mol o2 per mol for
//     glucose
 
//     MODEL_CONSTANT FloatType tau_1 = 1000.; // s
//     MODEL_CONSTANT FloatType tau_2 = 1000.; // s
 
//     MODEL_CONSTANT FloatType phi_max = (dl_max_ms * lin_density) * y_sx_1;
//     MODEL_CONSTANT FloatType phi_o2_max
//         = 10 * phi_max / MolarMassG * y_os_molar * MolarMassO2; // kgS/s
//     MODEL_CONSTANT FloatType nu_max_kg_s = dl_max_ms * lin_density;
 
//     MODEL_CONSTANT FloatType k_o
//         = 0.0001; // g/L: Anane et. al 2017 (Biochem. Eng. J) (g/g)
//     MODEL_CONSTANT FloatType k = 1e-2;
//     MODEL_CONSTANT FloatType k_perm = 1e-3;
//     MODEL_CONSTANT FloatType beta = 7;
//     MODEL_CONSTANT FloatType tau_ap_1 = 300;
//     MODEL_CONSTANT FloatType tau_ap_2 = 200;
//     MODEL_CONSTANT FloatType tau_ap_3 = 1000;
 
//     MODEL_CONSTANT auto length_c_dist
//         = MC::Distributions::TruncatedNormal<FloatType>(
//             l_c_m, l_c_m / 2., l_min_m, l_max_m); // use in out_str_l3
 
//     MODEL_CONSTANT FloatType adder_mean = 1.5e-6; // m
//     MODEL_CONSTANT auto adder_dist
//         = MC::Distributions::TruncatedNormal<FloatType>(
//             adder_mean,
//             adder_mean / 2.,
//             adder_mean / 20.,
//             adder_mean * 10); // use in out_str_l3
//     // Constants END
 
//     static MC::ContribIndexBounds get_bounds();
 
//     KOKKOS_INLINE_FUNCTION static double
//     mass(std::size_t idx, const SelfParticle& arr)
//     {
//       return GET_PROPERTY(Self::particle_var::length) * lin_density;
//     }
 
//     static std::vector<std::string_view>
//     names()
//     {
//       return { "age",  "length", "nu1",  "nu2",    "nu_eff_1", "nu_eff_2",
//                "a_p1", "a_p2",   "a_p3", "phi_o2", "phi_g",    "phi_pts" };
//     }
 
//     static std::vector<std::size_t>
//     get_number()
//     {
//       return { INDEX_FROM_ENUM(particle_var::age),
//                INDEX_FROM_ENUM(particle_var::length),
//                INDEX_FROM_ENUM(particle_var::nu1),
//                INDEX_FROM_ENUM(particle_var::nu2),
//                INDEX_FROM_ENUM(particle_var::nu_eff_1),
//                INDEX_FROM_ENUM(particle_var::nu_eff_2),
//                INDEX_FROM_ENUM(Uptakeparticle_var::ap_1),
//                INDEX_FROM_ENUM(Uptakeparticle_var::ap_2),
//                INDEX_FROM_ENUM(Uptakeparticle_var::ap_3),
//                INDEX_FROM_ENUM(particle_var::contrib_phi_o2),
//                INDEX_FROM_ENUM(particle_var::contrib_phi_s),
//                INDEX_FROM_ENUM(particle_var::phi_pts) };
//     }
 
//     static KOKKOS_INLINE_FUNCTION void
//     preinit()
//     {
//       Kokkos::printf("[Model]: PRENINIT:BEGIN\r\n");
//       //      Kokkos::printf("[Model]: phi_max:%.12f\r\n", phi_pts_max *
//       1e12);
//       //      Kokkos::printf("[Model]: PRENINIT:END\r\n");
//     }
 
//     KOKKOS_INLINE_FUNCTION static void init(const MC::pool_type& random_pool,
//                                             std::size_t idx,
//                                             const SelfParticle& arr);
 
//     KOKKOS_INLINE_FUNCTION static MC::Status
//     update(const MC::pool_type& random_pool,
//            FloatType d_t,
//            std::size_t idx,
//            const SelfParticle& arr,
//            const MC::LocalConcentration& c);
 
//     KOKKOS_INLINE_FUNCTION static void
//     division(const MC::pool_type& random_pool,
//              std::size_t idx,
//              std::size_t idx2,
//              const SelfParticle& arr,
//              const SelfParticle& buffer_arr);
 
//     KOKKOS_INLINE_FUNCTION static void
//     contribution(std::size_t idx,
//                  std::size_t position,
//                  double weight,
//                  const SelfParticle& arr,
//                  const MC::ContributionView& contributions);
//   };
 
//   CHECK_MODEL(TwoMetaNb)
 
//   KOKKOS_INLINE_FUNCTION void
//   TwoMetaNb::init([[maybe_unused]] const MC::pool_type& random_pool,
//                   std::size_t idx,
//                   const SelfParticle& arr)
//   {
 
//     // auto& v = init_uptake_cst;
//     constexpr auto local_ac = adder_dist;
 
//     constexpr auto length_dist =
//     MC::Distributions::TruncatedNormal<FloatType>(
//         l_c_m / 2, l_c_m / 5., l_min_m, l_max_m);
 
//     constexpr auto mu_nu_dist = nu_max_kg_s * 0.1;
//     constexpr auto nu_1_initial_dist
//         = MC::Distributions::TruncatedNormal<FloatType>(
//             mu_nu_dist, mu_nu_dist / 7., 0.,
//             static_cast<double>(nu_max_kg_s));
 
//     auto gen = random_pool.get_state();
//     auto l = length_dist.draw(gen);
 
//     GET_PROPERTY(Self::particle_var::length) = l;
//     GET_PROPERTY(Self::particle_var::l_cp) = l + local_ac.draw(gen);
//     GET_PROPERTY(particle_var::nu1) = nu_1_initial_dist.draw(gen);
//     random_pool.free_state(gen);
//     GET_PROPERTY(particle_var::contrib_phi_s) = 0;
 
//     Uptake<Self>::init(random_pool, idx, arr);
//   }
 
//   KOKKOS_INLINE_FUNCTION MC::Status
//   TwoMetaNb::update([[maybe_unused]] const MC::pool_type& random_pool,
//                     FloatType d_t,
//                     std::size_t idx,
//                     const SelfParticle& arr,
//                     const MC::LocalConcentration& concentrations)
//   {
//     const auto phi_s
//         = Uptake<Self>::uptake_step(phi_max,
//                                     d_t,
//                                     idx,
//                                     arr,
//                                     concentrations,
//                                     &GET_PROPERTY(Self::particle_var::phi_pts));
 
//     const auto o = Kokkos::max(static_cast<FloatType>(concentrations(1)),
//     0.F);
 
//     const auto phi_o2 = (phi_o2_max)*o / (o + k_o); // gO2/s
 
//     const auto nu_1_star
//         = y_sx_1 * MolarMassG
//           * Kokkos::min(phi_s / MolarMassG,
//                         phi_o2 / MolarMassO2 / y_os_molar); // gX/s
 
//     const auto s_1_star = (1.F / y_sx_1 * nu_1_star);
 
//     const auto phi_s_residual_1_star = Kokkos::max(phi_s - s_1_star, 0.F);
//     KOKKOS_ASSERT(phi_s_residual_1_star >= 0.F);
 
//     const auto nu_2_star = y_sx_2 * phi_s_residual_1_star; // gX/s
 
//     auto& nu_eff_1 = GET_PROPERTY(Self::particle_var::nu_eff_1);
//     auto& nu_eff_2 = GET_PROPERTY(Self::particle_var::nu_eff_2);
//     auto& nu_1 = GET_PROPERTY(Self::particle_var::nu1);
//     auto& nu_2 = GET_PROPERTY(Self::particle_var::nu2);
 
//     nu_eff_1 = Kokkos::min(nu_1_star, nu_1); // gX/s
 
//     const auto s_1 = (1.F / y_sx_1 * nu_eff_1);
//     const auto phi_s_residual_1 = Kokkos::max(phi_s - s_1, 0.F);
 
//     nu_eff_2 = Kokkos::min(y_sx_2 * phi_s_residual_1,
//                            nu_2); // gX/s
 
//     const auto sum_nu = (nu_eff_1 + nu_eff_2);
 
//     const auto s_growth = s_1 + (1 / y_sx_2 * nu_eff_2);
 
//     const auto s_overflow = phi_s - s_growth;
 
//     KOKKOS_ASSERT(nu_eff_1 >= 0.F);
//     KOKKOS_ASSERT(nu_eff_2 >= 0.F);
 
//     // CONTRIBS
//     GET_PROPERTY(Self::particle_var::contrib_phi_s) = -phi_s;
//     GET_PROPERTY(Self::particle_var::contrib_phi_o2)
//         = -1.F
//           * ((1.F / y_sx_1 / MolarMassG * y_os_molar * MolarMassO2 *
//           nu_eff_1)
//              + 0.F * nu_eff_2);
 
//     GET_PROPERTY(Self::particle_var::contrib_phi_ac)
//         = nu_eff_2 / y_sx_2 * y_sa
//           + (s_overflow > 0. ? y_sa * (s_overflow) : 0);
 
//     // ODE
//     nu_1 += d_t * ((nu_1_star - nu_1) / tau_1);
 
//     nu_2 += d_t * ((nu_2_star - nu_2) / tau_2);
 
//     GET_PROPERTY(Self::particle_var::length) += d_t * (sum_nu / lin_density);
 
//     GET_PROPERTY(Self::particle_var::age) += d_t;
 
//     return (GET_PROPERTY(Self::particle_var::length)
//             > GET_PROPERTY(Self::particle_var::l_cp))
//                ? MC::Status::Division
//                : MC::Status::Idle;
//   }
 
//   KOKKOS_INLINE_FUNCTION void
//   TwoMetaNb::division(const MC::pool_type& random_pool,
//                       std::size_t idx,
//                       std::size_t idx2,
//                       const SelfParticle& arr,
//                       const SelfParticle& child_buffer_arr)
//   {
//     constexpr FloatType half = 0.5;
//     constexpr auto local_ac = adder_dist;
//     const FloatType new_current_length
//         = GET_PROPERTY(particle_var::length) / static_cast<FloatType>(2.);
 
//     GET_PROPERTY(Self::particle_var::length) = new_current_length;
//     GET_PROPERTY(Self::particle_var::age) = 0;
 
//     GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::length)
//         = new_current_length;
 
//     GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::age) = 0;
 
//     const auto nu_1_o
//         = GET_PROPERTY_FROM(idx, arr, Self::particle_var::nu_eff_1);
//     const auto nu_2_o
//         = GET_PROPERTY_FROM(idx, arr, Self::particle_var::nu_eff_2);
 
//     auto gen = random_pool.get_state();
 
//     GET_PROPERTY(Self::particle_var::l_cp)
//         = new_current_length + local_ac.draw(gen);
 
//     if (nu_1_o != 0)
//     {
//       GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::nu1)
//           = MC::Distributions::TruncatedNormal<FloatType>::draw_from(
//               gen, nu_1_o, nu_1_o * half, 0.F, 1.F);
//     }
 
//     if (nu_2_o != 0)
//     {
//       GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::nu2)
//           = MC::Distributions::TruncatedNormal<FloatType>::draw_from(
//               gen, nu_2_o, nu_2_o * half, 0.F, 1.F);
//     }
 
//     GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::l_cp)
//         = new_current_length + local_ac.draw(gen);
//     random_pool.free_state(gen);
 
//     Uptake<Self>::division(random_pool, idx, idx2, arr, child_buffer_arr);
//   }
 
//   inline MC::ContribIndexBounds
//   TwoMetaNb::get_bounds()
//   {
//     const int begin = INDEX_FROM_ENUM(Self::particle_var::contrib_phi_s);
//     const int end = INDEX_FROM_ENUM(Self::particle_var::contrib_phi_ac);
//     assert(begin < end);
//     return { .begin = begin, .end = end };
//   }
//   // KOKKOS_INLINE_FUNCTION void
//   // TwoMetaNb::contribution([[maybe_unused]] std::size_t idx,
//   //                         std::size_t position,
//   //                         double weight,
//   //                         [[maybe_unused]] const SelfParticle& arr,
//   //                         const MC::ContributionView& contributions)
//   // {
//   //   auto access = contributions.access();
//   //   access(0, position) +=
//   //       weight * GET_PROPERTY(Self::particle_var::contrib_phi_s); //
//   NOLINT
//   //   access(1, position) +=
//   //       weight * GET_PROPERTY(Self::particle_var::contrib_phi_o2); //
//   NOLINT
//   //   access(2, position) +=
//   //       weight * GET_PROPERTY(Self::particle_var::contrib_phi_ac); //
//   NOLINT
//   // }
 
//   static_assert(HasExportProperties<TwoMetaNb>, "ee");
 
// } // namespace Models
 
// #endif //__TWOMETA_MODEL_NB_HPP__