rscm 0.2.3

Rust Simple Climate Model - A framework for building reduced-complexity climate models
Documentation 
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

#![allow(dead_code)]

use ode_solvers::*;
use std::collections::HashMap;
use std::sync::Arc;

use rscm_core::component::{
    Component, InputState, OutputState, RequirementDefinition, RequirementType,
};
use rscm_core::errors::RSCMResult;
use rscm_core::ivp::{IVPBuilder, IVP};
use rscm_core::timeseries::{FloatValue, Time};
use serde::{Deserialize, Serialize};

// Define some types that are used by OdeSolvers
type ModelState = Vector3<FloatValue>;

#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct TwoLayerComponentParameters {
    pub lambda0: FloatValue,
    pub a: FloatValue,
    pub efficacy: FloatValue,
    pub eta: FloatValue,
    pub heat_capacity_surface: FloatValue,
    pub heat_capacity_deep: FloatValue,
}

#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct TwoLayerComponent {
    parameters: TwoLayerComponentParameters,
}

// Create the set of ODEs to represent the two layer model
impl IVP<Time, ModelState> for TwoLayerComponent {
    fn calculate_dy_dt(
        &self,
        _t: Time,
        input_state: &InputState,
        y: &ModelState,
        dy_dt: &mut ModelState,
    ) {
        let temperature_surface = y[0];
        let temperature_deep = y[1];
        let erf = input_state.get_latest("Effective Radiative Forcing");

        let temperature_difference = temperature_surface - temperature_deep;

        let lambda_eff = self.parameters.lambda0 - self.parameters.a * temperature_surface;
        let heat_exchange_surface =
            self.parameters.efficacy * self.parameters.eta * temperature_difference;
        let dtemperature_surface_dt =
            (erf - lambda_eff * temperature_surface - heat_exchange_surface)
                / self.parameters.heat_capacity_surface;

        let heat_exchange_deep = self.parameters.eta * temperature_difference;
        let dtemperature_deep_dt = heat_exchange_deep / self.parameters.heat_capacity_deep;

        dy_dt[0] = dtemperature_surface_dt;
        dy_dt[1] = dtemperature_deep_dt;
        dy_dt[2] = self.parameters.heat_capacity_surface * dtemperature_surface_dt
            + self.parameters.heat_capacity_deep * dtemperature_deep_dt;
    }
}

impl TwoLayerComponent {
    pub fn from_parameters(parameters: TwoLayerComponentParameters) -> Self {
        Self { parameters }
    }
}

#[typetag::serde]
impl Component for TwoLayerComponent {
    fn definitions(&self) -> Vec<RequirementDefinition> {
        vec![
            RequirementDefinition::new(
                "Effective Radiative Forcing",
                "W/m^2",
                RequirementType::Input,
            ),
            RequirementDefinition::new("Surface Temperature", "K", RequirementType::Output),
        ]
    }

    fn solve(
        &self,
        t_current: Time,
        t_next: Time,
        input_state: &InputState,
    ) -> RSCMResult<OutputState> {
        let erf = input_state.get_latest("Effective Radiative Forcing");

        let y0 = ModelState::new(0.0, 0.0, 0.0);

        let solver = IVPBuilder::new(Arc::new(self.to_owned()), input_state, y0);
        println!("Solving {:?} with state: {:?}", self, input_state);

        let mut solver = solver.to_rk4(t_current, t_next, 0.1);
        let stats = solver.integrate().expect("Failed solving");

        let results = solver.results();

        println!("Stats {:?}", stats);
        println!("Results {:?}", results);

        // Create the solver

        Ok(HashMap::from([(
            "Surface Temperature".to_string(),
            erf * self.parameters.lambda0,
        )]))
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use numpy::array;
    use rscm_core::model::extract_state;
    use rscm_core::timeseries::Timeseries;
    use rscm_core::timeseries_collection::{TimeseriesCollection, VariableType};

    #[test]
    fn it_works() {
        // Solve the two layer component in isolation
        let component = TwoLayerComponent::from_parameters(TwoLayerComponentParameters {
            lambda0: 0.5,
            a: 0.01,
            efficacy: 0.5,
            eta: 0.1,
            heat_capacity_surface: 1.0,
            heat_capacity_deep: 100.0,
        });

        let mut ts_collection = TimeseriesCollection::new();
        ts_collection.add_timeseries(
            "Effective Radiative Forcing".to_string(),
            Timeseries::from_values(
                array![1.0, 1.5, 2.0, 2.0],
                array![1848.0, 1849.0, 1850.0, 1900.0],
            ),
            VariableType::Exogenous,
        );

        let input_state = extract_state(&ts_collection, component.input_names(), 1848.0);

        // Create the solver
        let output_state = component.solve(1848.0, 1849.0, &input_state);

        println!("Output: {:?}", output_state);
        let output_state = output_state.unwrap();
        assert_eq!(*output_state.get("Surface Temperature").unwrap(), 0.5);
    }
}