1+ #include < algorithm>
2+ #include < cma_utils/alias.hpp>
3+ #include < limits>
4+ namespace CmaUtils
5+ {
6+
7+ // double get_min_residence_time(const TransitionnerPtrType& iterator)
8+ // {
9+ // const std::size_t n_states = iterator->size();
10+ // double min_residence_time = std::numeric_limits<double>::max();
11+ // for (std::size_t i_state = 0; i_state < n_states; ++i_state)
12+ // {
13+ // // Naive min-research, as COO is not sorted we can´t do better
14+ // without
15+ // // convert COO to CSR/CSC or sort
16+ // // + function exececuted only once
17+ // const auto state = iterator->get_at(i_state);
18+ // const auto liquid = state->get_liquid();
19+ // const auto coo_matrix = liquid->transition();
20+ // const auto liquid_volumes = liquid->volume();
21+ // const auto data = coo_matrix->values();
22+ // const auto rows = coo_matrix->row_indices();
23+ // const auto cols = coo_matrix->col_indices();
24+
25+ // for (std::size_t k = 0; k < data.size(); ++k)
26+ // {
27+ // // Only way to find diagonal with not sorted COO
28+ // if (rows[k] == cols[k])
29+ // {
30+ // const auto volume = liquid_volumes[rows[k]];
31+ // if (volume != 0.)
32+ // {
33+ // double residence_time =
34+ // -data[k] / volume; // minus cause transition has negative
35+ // // diagonal (leaving flow)
36+ // min_residence_time = std::min(residence_time,
37+ // min_residence_time);
38+ // }
39+ // }
40+ // }
41+ // }
42+ // return min_residence_time;
43+ // }
44+
45+ // Does this function should be moved into RCMTool ?
46+ double get_min_residence_time (const TransitionnerPtrType& iterator) noexcept
47+ {
48+ /*
49+ This new implementation only determine minimum residence time defined as
50+ min(tau)=min(qi/vi) where qi=sum(flow) in compartment i First impl
51+ determined the smaller flow min(flowi/vi) which leads to smaller value
52+ */
53+ const std::size_t n_states = iterator->size ();
54+ double min_residence_time = std::numeric_limits<double >::max ();
55+ for (std::size_t i_state = 0 ; i_state < n_states; ++i_state)
56+ {
57+ const auto state = iterator->get_at (i_state);
58+ const auto liquid = state->get_liquid ();
59+ const auto out_flows = liquid->out_flows ();
60+ const auto liquid_volumes = liquid->volume ();
61+ KOKKOS_ASSERT (liquid_volumes.size () == out_flows.size ());
62+ for (std::size_t k = 0 ; k < out_flows.size (); ++k)
63+ {
64+ const auto volume = liquid_volumes[k];
65+ if (volume > 0 .)
66+ {
67+ const double residence_time = out_flows[k] / volume;
68+ min_residence_time = std::min (residence_time, min_residence_time);
69+ }
70+ }
71+ }
72+ return min_residence_time;
73+ }
74+
75+ } // namespace CmaUtils
0 commit comments