1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50

//! Reactor structure.

use crate::chem::Rate;
use ndarray::{Array1, Array2};

/// Total reactor structure.
#[derive(Debug)]
pub struct Reactor {
    /// Rates.
    rates: Array1<Rate>,
    /// Stoichiometric coefficents.
    coeffs: Array2<f64>,
}

impl Reactor {
    /// Construct a new instance.
    #[inline]
    #[must_use]
    pub fn new(rates: Array1<Rate>, coeffs: Array2<f64>) -> Self {
        debug_assert!(!rates.is_empty());
        debug_assert!(rates.len() == coeffs.nrows());

        Self { rates, coeffs }
    }

    /// Calculate the reaction rates.
    #[inline]
    #[must_use]
    fn rates(&self, concs: &Array1<f64>) -> Array1<f64> {
        debug_assert!(concs.len() == self.coeffs.ncols());

        self.rates.iter().map(|r| r.rate(concs)).collect()
    }

    /// Calculate the overall rate of change given the current concentrations.
    #[inline]
    #[must_use]
    pub fn deltas(&self, concs: &Array1<f64>) -> Array1<f64> {
        debug_assert!(concs.len() == self.coeffs.ncols());

        let rates = self.rates(concs);

        let mut deltas = Array1::zeros(concs.len());
        for (coeffs, rate) in self.coeffs.outer_iter().zip(&rates) {
            deltas += &(&coeffs * *rate);
        }

        deltas
    }
}