Commit f0e76a99 authored by Turner, Sean's avatar Turner, Sean
Browse files

hard-coded single reservoir

parent a495214b
Loading
Loading
Loading
Loading
+42 −3
Original line number Diff line number Diff line
@@ -3,14 +3,53 @@ use pyo3::exceptions::PyValueError;

/// Simulates a cascade of reservoirs
#[pyfunction]
fn simulate_cascade(inflow: Vec<f64>, release: Vec<f64>) -> PyResult<Vec<f64>> {
fn simulate_cascade(inflow: Vec<f64>, release: Vec<f64>) -> PyResult<(Vec<f64>, Vec<f64>, Vec<f64>)> {

    // Reservoir specifications
    const CAPACITY: f64 = 100.0;  // Maximum storage in MCM (Million cubic meters)
    const INITIAL_STORAGE: f64 = 50.0;  // Initial storage in MCM
    const MAX_RELEASE: f64 = 15.0;  // Maximum release per timestep in MCM
    const MIN_RELEASE: f64 = 3.0;  // Minimum release per timestep in MCM

    if inflow.len() != release.len() {
        return Err(PyValueError::new_err("Input vectors must have the same length"));
    }

    // placeholder: sums inflow and release to return single vector
    let n = inflow.len();
    let mut storage = vec![0.0; n];  // Storage at end of timestep
    let mut actual_release = vec![0.0; n];  // Actual release for timestep
    let mut spill = vec![0.0; n];  // Spill for timestep

    // Initialize first timestep storage
    storage[0] = INITIAL_STORAGE;

    for t in 0..n {
        let available_water = if t == 0 {
            INITIAL_STORAGE + inflow[t]
        } else {
            storage[t-1] + inflow[t]
        };

        // Calculate actual release (constrained by available water and max/min release)
        actual_release[t] = release[t]
            .max(MIN_RELEASE)  // Ensure we meet minimum release
            .min(MAX_RELEASE)  // Cannot exceed maximum release
            .min(available_water);  // Cannot release more than available

        // Calculate new storage
        let potential_storage = available_water - actual_release[t];

        // Calculate spill if storage would exceed capacity
        if potential_storage > CAPACITY {
            storage[t] = CAPACITY;
            spill[t] = potential_storage - CAPACITY;
        } else {
            storage[t] = potential_storage;
            spill[t] = 0.0;
        }
    }

    Ok(inflow.iter().zip(release.iter()).map(|(x, y)| x + y).collect())
    Ok((storage, actual_release, spill))

}