BioCMAMC-ST
two_meta.hpp
1#ifndef __TWOMETA_MODEL_HPP__
2#define __TWOMETA_MODEL_HPP__
3
4#include "Kokkos_Printf.hpp"
5#include "common/traits.hpp"
6#include "mc/macros.hpp"
7#include "models/utils.hpp"
8#include <mc/prng/prng_extension.hpp>
9#include <mc/traits.hpp>
10#include <models/uptake.hpp>
11#include <string_view>
12
13namespace
14{
15 // template <FloatingPointType F> static F consteval get_phi_s_max(F density, F dl)
16 // {
17 // // dl and density must be same unit, dl*density -> mass and y is mass yield
18 // return (dl * density) * 0.5;
19 // }
20} // namespace
21
22namespace Models
23{
24
25 struct TwoMeta
26 {
27 using uniform_weight = std::true_type; // Using type alias
28 using Self = TwoMeta;
29 using FloatType = float;
30 using Config = std::nullopt_t;
31
32 enum class particle_var : int
33 {
34 age = INDEX_FROM_ENUM(Uptakeparticle_var::COUNT),
35 length,
36 nu1,
37 nu2,
38 l_cp,
39 nu_eff_1, // This is not itself a model property but stored to be exported
40 nu_eff_2, // This is not itself a model property but stored to be exported
41 // TODO FIND BETTER WAY TO STORE/GET CONTRIBUTIONS
45 COUNT
46 };
47
48 static constexpr std::size_t n_var = INDEX_FROM_ENUM(particle_var::COUNT);
49 static constexpr std::string_view name = "simple";
51
52 MODEL_CONSTANT FloatType l_max_m = 5e-6; // m
53 MODEL_CONSTANT FloatType l_c_m = 3e-6; // m
54 MODEL_CONSTANT FloatType d_m = 0.6e-6; // m
55 MODEL_CONSTANT FloatType l_min_m = 0.9e-6; // m
56 MODEL_CONSTANT FloatType lin_density = c_linear_density(static_cast<FloatType>(1000), d_m);
59
60 MODEL_CONSTANT FloatType y_sx_1 = 1. / 2.217737e+00; // Mode 1 S to X yield (mass)
61 MODEL_CONSTANT FloatType y_sx_2 = y_sx_1 / 3.; // Mode 2 S to X yield (mass)
62 MODEL_CONSTANT FloatType y_sa = 0.8; // S to A yield (mass)
63 MODEL_CONSTANT FloatType y_os_molar = 3; // 3 mol o2 per mol for glucose
64 MODEL_CONSTANT FloatType k_o = 0.0001; // g/L: Anane et. al 2017 (Biochem. Eng. J) (g/g)
65 MODEL_CONSTANT FloatType dl_max_ms = 8 * 2e-10; // m/s https://doi.org/10.7554/eLife.67495;
66 MODEL_CONSTANT FloatType tau_1 = 1000.; // s
67 MODEL_CONSTANT FloatType tau_2 = 1000.; // s
68
70 MODEL_CONSTANT FloatType phi_perm_max = phi_s_max / 40.; // kgS/
71 MODEL_CONSTANT FloatType phi_o2_max =
72 10 * phi_s_max / MolarMassG * y_os_molar * MolarMassO2; // kgS/s
73 MODEL_CONSTANT float nu_max_kg_s = dl_max_ms * lin_density;
74
75 // MODEL_CONSTANT auto length_c_dist =
76 // MC::Distributions::TruncatedNormal<FloatType>(l_c_m, l_c_m / 7., l_min_m, l_max_m); use
77 // in out_str_l2
78
80 l_c_m, l_c_m / 2., l_min_m, l_max_m); // use in out_str_l3
81
82 // MODEL_CONSTANT auto length_c_dist =
83 // MC::Distributions::TruncatedNormal<FloatType>(1.5*l_c_m, l_c_m / 7., 3*l_min_m, l_max_m);
84
85 KOKKOS_INLINE_FUNCTION static void
86 init(const MC::KPRNG::pool_type& random_pool, std::size_t idx, const SelfParticle& arr);
87
88 KOKKOS_INLINE_FUNCTION static MC::Status update(const MC::KPRNG::pool_type& random_pool,
89 FloatType d_t,
90 std::size_t idx,
91 const SelfParticle& arr,
92 const MC::LocalConcentration& c);
93
94 KOKKOS_INLINE_FUNCTION static void division(const MC::KPRNG::pool_type& random_pool,
95 std::size_t idx,
96 std::size_t idx2,
97 const SelfParticle& arr,
98 const SelfParticle& buffer_arr);
99
100 KOKKOS_INLINE_FUNCTION static void contribution(std::size_t idx,
101 std::size_t position,
102 double weight,
103 const SelfParticle& arr,
104 const MC::ContributionView& contributions);
105
106 KOKKOS_INLINE_FUNCTION static double mass(std::size_t idx, const SelfParticle& arr)
107 {
108 return GET_PROPERTY(Self::particle_var::length) * lin_density;
109 }
110
111 static std::vector<std::string_view> names()
112 {
113 return {"age",
114 "length",
115 "nu1",
116 "nu2",
117 "nu_eff_1",
118 "nu_eff_2",
119 "a_perm",
120 "a_pts",
121 "phi_o2",
122 "phi_g"};
123 }
124
125 static std::vector<std::size_t> get_number()
126 {
127 return {INDEX_FROM_ENUM(particle_var::age),
128 INDEX_FROM_ENUM(particle_var::length),
129 INDEX_FROM_ENUM(particle_var::nu1),
130 INDEX_FROM_ENUM(particle_var::nu2),
131 INDEX_FROM_ENUM(particle_var::nu_eff_1),
132 INDEX_FROM_ENUM(particle_var::nu_eff_2),
133 INDEX_FROM_ENUM(Uptakeparticle_var::a_permease),
134 INDEX_FROM_ENUM(Uptakeparticle_var::a_pts),
135 INDEX_FROM_ENUM(particle_var::contrib_phi_o2),
136 INDEX_FROM_ENUM(particle_var::contrib_phi_s)};
137 }
138 };
139
140 CHECK_MODEL(TwoMeta)
141
142 KOKKOS_INLINE_FUNCTION void
143 TwoMeta::init([[maybe_unused]] const MC::KPRNG::pool_type& random_pool,
144 std::size_t idx,
145 const SelfParticle& arr)
146 {
147 constexpr auto local_lc = length_c_dist;
148 constexpr auto length_dist =
150
151 constexpr auto mu_nu_dist = nu_max_kg_s * 0.1;
152 constexpr auto nu_1_initial_dist = MC::Distributions::TruncatedNormal<float>(
153 mu_nu_dist, mu_nu_dist / 7., 0., static_cast<double>(nu_max_kg_s));
154
155 auto gen = random_pool.get_state();
156 GET_PROPERTY(Self::particle_var::length) = length_dist.draw(gen);
157 GET_PROPERTY(Self::particle_var::l_cp) = local_lc.draw(gen);
158 GET_PROPERTY(particle_var::nu1) = nu_1_initial_dist.draw(gen);
159 random_pool.free_state(gen);
160 GET_PROPERTY(particle_var::contrib_phi_s) = 0;
161 Uptake<Self>::init(random_pool, idx, arr);
162 }
163
164 KOKKOS_INLINE_FUNCTION MC::Status TwoMeta::update(const MC::KPRNG::pool_type& random_pool,
165 FloatType d_t,
166 std::size_t idx,
167 const SelfParticle& arr,
168 const MC::LocalConcentration& concentrations)
169 {
170 (void)random_pool;
171 const auto phi_s =
172 Uptake<Self>::uptake_step(phi_s_max, phi_perm_max, d_t, idx, arr, concentrations);
173
174 const auto o = Kokkos::max(static_cast<float>(concentrations(1)), 0.F);
175
176 const float phi_o2 = (phi_o2_max)*o / (o + k_o); // gO2/s
177
178 const float nu_1_star = y_sx_1 * MolarMassG *
179 Kokkos::min(phi_s / MolarMassG,
180 phi_o2 / MolarMassO2 / y_os_molar); // gX/s
181
182 const float s_1_star = (1 / y_sx_1 * nu_1_star);
183
184 const float phi_s_residual_1_star = Kokkos::max(phi_s - s_1_star, 0.F);
185 KOKKOS_ASSERT(phi_s_residual_1_star >= 0.F);
186
187 const float nu_2_star = y_sx_2 * phi_s_residual_1_star; // gX/s
188
189 GET_PROPERTY(Self::particle_var::nu_eff_1) =
190 Kokkos::min(nu_1_star, GET_PROPERTY(Self::particle_var::nu1)); // gX/s
191
192 const float s_1 = (1 / y_sx_1 * GET_PROPERTY(Self::particle_var::nu_eff_1));
193 const float phi_s_residual_1 = Kokkos::max(phi_s - s_1, 0.F);
194 GET_PROPERTY(Self::particle_var::nu_eff_2) =
195 Kokkos::min(y_sx_2 * phi_s_residual_1, GET_PROPERTY(Self::particle_var::nu2)); // gX/s
196
197 const float s_growth = s_1 + (1 / y_sx_2 * GET_PROPERTY(Self::particle_var::nu_eff_2));
198
199 const float s_overflow = phi_s - s_growth;
200
201 KOKKOS_ASSERT(GET_PROPERTY(Self::particle_var::nu_eff_1) >= 0.F);
202 KOKKOS_ASSERT(GET_PROPERTY(Self::particle_var::nu_eff_2) >= 0.F);
203
204 // CONTRIBS
205 GET_PROPERTY(Self::particle_var::contrib_phi_s) = -phi_s;
207 -1 * ((1. / y_sx_1 / MolarMassG * y_os_molar * MolarMassO2 *
208 GET_PROPERTY(Self::particle_var::nu_eff_1)) +
209 0. * GET_PROPERTY(Self::particle_var::nu_eff_2));
210
212 GET_PROPERTY(Self::particle_var::nu_eff_2) / y_sx_2 * y_sa +
213 (s_overflow > 0. ? y_sa * (s_overflow) : 0);
214
215 // ODE
216 GET_PROPERTY(Self::particle_var::nu1) +=
217 static_cast<float>(d_t) * ((nu_1_star - GET_PROPERTY(Self::particle_var::nu1)) / tau_1);
218
219 GET_PROPERTY(Self::particle_var::nu2) +=
220 static_cast<float>(d_t) * ((nu_2_star - GET_PROPERTY(Self::particle_var::nu2)) / tau_2);
221
222 const auto sum_nu =
223 (GET_PROPERTY(Self::particle_var::nu_eff_1) + GET_PROPERTY(Self::particle_var::nu_eff_2));
224
225 GET_PROPERTY(Self::particle_var::length) += static_cast<float>(d_t) * (sum_nu / lin_density);
226
227 GET_PROPERTY(Self::particle_var::age) += d_t;
228
229 return (GET_PROPERTY(Self::particle_var::length) > GET_PROPERTY(Self::particle_var::l_cp))
232 }
233
234 KOKKOS_INLINE_FUNCTION void TwoMeta::division(const MC::KPRNG::pool_type& random_pool,
235 std::size_t idx,
236 std::size_t idx2,
237 const SelfParticle& arr,
238 const SelfParticle& child_buffer_arr)
239 {
240 constexpr auto local_lc = length_c_dist;
241 const FloatType new_current_length =
242 GET_PROPERTY(particle_var::length) / static_cast<FloatType>(2.);
243
244 GET_PROPERTY(Self::particle_var::length) = new_current_length;
245 GET_PROPERTY(Self::particle_var::age) = 0;
246
247 GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::length) = new_current_length;
248 GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::age) = 0;
249
250 auto nu_1_o = GET_PROPERTY_FROM(idx, arr, Self::particle_var::nu_eff_1);
251 auto nu_2_o = GET_PROPERTY_FROM(idx, arr, Self::particle_var::nu_eff_2);
252 auto gen = random_pool.get_state();
253 GET_PROPERTY(Self::particle_var::l_cp) = local_lc.draw(gen);
254 if (nu_1_o != 0)
255 {
256 GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::nu1) =
257 MC::Distributions::TruncatedNormal<FloatType>::draw_from(gen, nu_1_o, nu_1_o / 2., 0, 1.);
258 }
259 if (nu_2_o != 0)
260 {
261 GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::nu2) =
262 MC::Distributions::TruncatedNormal<FloatType>::draw_from(gen, nu_2_o, nu_2_o / 2., 0, 1.);
263 }
264
265 GET_PROPERTY_FROM(idx2, child_buffer_arr, Self::particle_var::l_cp) = local_lc.draw(gen);
266 random_pool.free_state(gen);
267
268 Uptake<Self>::division(random_pool, idx, idx2, arr, child_buffer_arr);
269 }
270
271 KOKKOS_INLINE_FUNCTION void TwoMeta::contribution([[maybe_unused]] std::size_t idx,
272 std::size_t position,
273 double weight,
274 [[maybe_unused]] const SelfParticle& arr,
275 const MC::ContributionView& contributions)
276 {
277 auto access = contributions.access();
278 access(0, position) += weight * GET_PROPERTY(Self::particle_var::contrib_phi_s); // NOLINT
279 access(1, position) += weight * GET_PROPERTY(Self::particle_var::contrib_phi_o2); // NOLINT
280 access(2, position) += weight * GET_PROPERTY(Self::particle_var::contrib_phi_ac); // NOLINT
281 }
282
283 static_assert(HasExportProperties<TwoMeta>, "ee");
284
285} // namespace Models
286
287#endif
Kokkos::Random_XorShift1024_Pool< Kokkos::DefaultExecutionSpace > pool_type
Definition prng.hpp:17
Model that can export properties.
Definition traits.hpp:159
Kokkos::Subview< KernelConcentrationType, int, decltype(Kokkos::ALL)> LocalConcentration
Definition alias.hpp:42
Kokkos::Experimental::ScatterView< double **, Kokkos::LayoutRight > ContributionView
Definition alias.hpp:31
Kokkos::View< F *[Nd]> ParticlesModel
Definition traits.hpp:34
Status
Definition alias.hpp:11
constexpr T glucose
Definition utils.hpp:21
constexpr T dioxygen
Definition utils.hpp:22
Models definition.
Definition config_loader.hpp:8
@ a_permease
Definition uptake.hpp:15
@ a_pts
Definition uptake.hpp:14
@ COUNT
Definition uptake.hpp:17
KOKKOS_INLINE_FUNCTION consteval F c_linear_density(F rho, F d)
Definition utils.hpp:39
Represents a TruncatedNormal (Gaussian) probability distribution.
Definition prng_extension.hpp:263
static KOKKOS_INLINE_FUNCTION F draw_from(Kokkos::Random_XorShift1024< DeviceType > &gen, F mu, F sigma, F lower, F upper)
Definition prng_extension.hpp:287
Definition two_meta.hpp:26
MODEL_CONSTANT FloatType l_c_m
Definition two_meta.hpp:53
MODEL_CONSTANT FloatType lin_density
Definition two_meta.hpp:56
static std::vector< std::size_t > get_number()
Definition two_meta.hpp:125
MODEL_CONSTANT FloatType phi_o2_max
Definition two_meta.hpp:71
MODEL_CONSTANT FloatType MolarMassG
Definition two_meta.hpp:57
static KOKKOS_INLINE_FUNCTION void contribution(std::size_t idx, std::size_t position, double weight, const SelfParticle &arr, const MC::ContributionView &contributions)
Definition two_meta.hpp:271
static constexpr std::size_t n_var
Definition two_meta.hpp:48
static constexpr std::string_view name
Definition two_meta.hpp:49
MODEL_CONSTANT FloatType l_max_m
Definition two_meta.hpp:52
MODEL_CONSTANT FloatType phi_perm_max
Definition two_meta.hpp:70
MODEL_CONSTANT FloatType k_o
Definition two_meta.hpp:64
MODEL_CONSTANT auto length_c_dist
Definition two_meta.hpp:79
static KOKKOS_INLINE_FUNCTION void init(const MC::KPRNG::pool_type &random_pool, std::size_t idx, const SelfParticle &arr)
Definition two_meta.hpp:143
MODEL_CONSTANT FloatType tau_2
Definition two_meta.hpp:67
float FloatType
Definition two_meta.hpp:29
MODEL_CONSTANT FloatType MolarMassO2
Definition two_meta.hpp:58
MODEL_CONSTANT FloatType y_sx_1
Definition two_meta.hpp:60
MODEL_CONSTANT FloatType l_min_m
Definition two_meta.hpp:55
static KOKKOS_INLINE_FUNCTION double mass(std::size_t idx, const SelfParticle &arr)
Definition two_meta.hpp:106
std::nullopt_t Config
Definition two_meta.hpp:30
MODEL_CONSTANT FloatType y_sa
Definition two_meta.hpp:62
MODEL_CONSTANT FloatType d_m
Definition two_meta.hpp:54
std::true_type uniform_weight
Definition two_meta.hpp:27
particle_var
Definition two_meta.hpp:33
MODEL_CONSTANT FloatType phi_s_max
Definition two_meta.hpp:69
MODEL_CONSTANT FloatType y_sx_2
Definition two_meta.hpp:61
static std::vector< std::string_view > names()
Definition two_meta.hpp:111
static KOKKOS_INLINE_FUNCTION MC::Status update(const MC::KPRNG::pool_type &random_pool, FloatType d_t, std::size_t idx, const SelfParticle &arr, const MC::LocalConcentration &c)
Definition two_meta.hpp:164
MC::ParticlesModel< Self::n_var, Self::FloatType > SelfParticle
Definition two_meta.hpp:50
MODEL_CONSTANT FloatType tau_1
Definition two_meta.hpp:66
MODEL_CONSTANT FloatType y_os_molar
Definition two_meta.hpp:63
MODEL_CONSTANT float nu_max_kg_s
Definition two_meta.hpp:73
MODEL_CONSTANT FloatType dl_max_ms
Definition two_meta.hpp:65
static KOKKOS_INLINE_FUNCTION void division(const MC::KPRNG::pool_type &random_pool, std::size_t idx, std::size_t idx2, const SelfParticle &arr, const SelfParticle &buffer_arr)
Definition two_meta.hpp:234
static KOKKOS_INLINE_FUNCTION void init(const MC::KPRNG::pool_type &random_pool, std::size_t idx, const SelfParticle &arr)
Definition uptake.hpp:68
static KOKKOS_INLINE_FUNCTION FloatType uptake_step(FloatType phi_pts_max, FloatType phi_permease_specific, FloatType d_t, std::size_t idx, const SelfParticle &arr, const MC::LocalConcentration &c)
Definition uptake.hpp:105
static KOKKOS_INLINE_FUNCTION void division(const MC::KPRNG::pool_type &random_pool, std::size_t idx, std::size_t idx2, const SelfParticle &arr, const SelfParticle &buffer_arr)
Definition uptake.hpp:174