Skip to content

Commit f546181

Browse files
LevanBokeriaLevanBokeria
authored
Merge pull request #70 from alan-turing-institute/66-improved-contraction-pressure
66 improved contraction pressure
2 parents 81fe830 + e9e0e26 commit f546181

File tree

7 files changed

+875
-3
lines changed

7 files changed

+875
-3
lines changed

‎sandbox/notebooks/generate_KMMPP_KGO.ipynb

Lines changed: 150 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,150 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Description\n",
8+
"\n",
9+
"This notebook is to generate known good outputs (KGOs) for solving KorakianitisMixedModelPP (KMMPP) after we have added a cycle step size argument to the `advance_cycle()` function of the solver.\n",
10+
"\n",
11+
"With different step sizes, there are slightly different solutions reached, so tests need these different KGOs to have as comparisons."
12+
]
13+
},
14+
{
15+
"cell_type": "code",
16+
"execution_count": 1,
17+
"metadata": {},
18+
"outputs": [],
19+
"source": [
20+
"import numpy as np\n",
21+
"import json\n",
22+
"\n",
23+
"from ModularCirc.Solver import Solver\n",
24+
"\n",
25+
"from ModularCirc.Models.KorakianitisMixedModelPP import KorakianitisMixedModelPP\n",
26+
"from ModularCirc.Models.KorakianitisMixedModel_parameters import KorakianitisMixedModel_parameters"
27+
]
28+
},
29+
{
30+
"cell_type": "code",
31+
"execution_count": null,
32+
"metadata": {},
33+
"outputs": [
34+
{
35+
"name": "stdout",
36+
"output_type": "stream",
37+
"text": [
38+
"Running simulation with cycle step size: 1\n",
39+
"Running simulation with cycle step size: 1\n",
40+
"Running simulation with cycle step size: 3\n",
41+
"Running simulation with cycle step size: 3\n",
42+
"Running simulation with cycle step size: 5\n",
43+
"Running simulation with cycle step size: 5\n",
44+
"Running simulation with cycle step size: 7\n",
45+
"Running simulation with cycle step size: 7\n"
46+
]
47+
}
48+
],
49+
"source": [
50+
"time_setup_dict = {\n",
51+
" 'name': 'TimeTest',\n",
52+
" 'ncycles': 40,\n",
53+
" 'tcycle': 1.0,\n",
54+
" 'dt': 0.001,\n",
55+
" 'export_min': 1\n",
56+
" }\n",
57+
"\n",
58+
"# Start the loop with different cycle step size to generate JSON files for each loop:\n",
59+
"cycle_step_sizes = [1, 3, 5, 7]\n",
60+
"\n",
61+
"# Dictionary to store the results for all cycle step sizes\n",
62+
"all_results = {\n",
63+
" \"metadata\": {\n",
64+
" \"description\": \"Results for different cycle_step_size values\",\n",
65+
" \"cycle_step_sizes\": cycle_step_sizes\n",
66+
" },\n",
67+
" \"results\": {}\n",
68+
"}\n",
69+
"\n",
70+
"# Loop through each cycle step size\n",
71+
"for i_cycle_step_size in cycle_step_sizes:\n",
72+
"\n",
73+
" print(f\"Running simulation with cycle step size: {i_cycle_step_size}\")\n",
74+
"\n",
75+
" # Initializing the parameter object\n",
76+
" parobj = KorakianitisMixedModel_parameters()\n",
77+
"\n",
78+
" # Initializing the model \n",
79+
" model = KorakianitisMixedModelPP(time_setup_dict=time_setup_dict, \n",
80+
" parobj=parobj, \n",
81+
" suppress_printing=True)\n",
82+
"\n",
83+
" # Initializing the solver\n",
84+
" solver = Solver(model=model)\n",
85+
"\n",
86+
" # Initializing the parameter object\n",
87+
" parobj = KorakianitisMixedModel_parameters()\n",
88+
"\n",
89+
" # Initializing the model \n",
90+
" model = KorakianitisMixedModelPP(time_setup_dict=time_setup_dict, \n",
91+
" parobj=parobj, \n",
92+
" suppress_printing=True)\n",
93+
"\n",
94+
" # Initializing the solver\n",
95+
" solver = Solver(model=model)\n",
96+
"\n",
97+
" # Solver is being setup: switching off console printing and setting the solver method to \"LSODA\"\n",
98+
" solver.setup(suppress_output=True, \n",
99+
" method='LSODA',\n",
100+
" step=i_cycle_step_size)\n",
101+
"\n",
102+
" # Running the model\n",
103+
" solver.solve()\n",
104+
"\n",
105+
" # Define the indexes of the equivalent to the last cycles\n",
106+
" tind_fin = np.arange(start=model.time_object.n_t-model.time_object.n_c,\n",
107+
" stop=model.time_object.n_t)\n",
108+
"\n",
109+
" # Dictionary to store the final cycle values for each component for the current step size\n",
110+
" final_cycle_values = {}\n",
111+
"\n",
112+
" # From each of the components, retrieve the volume (<V>), pressure (<P_i>) and flow (<Q_i>)\n",
113+
" for key, value in model.components.items():\n",
114+
" final_cycle_values[key] = {\n",
115+
" 'V': value.V.values[tind_fin].mean(),\n",
116+
" 'P_i': value.P_i.values[tind_fin].mean(),\n",
117+
" 'Q_i': value.Q_i.values[tind_fin].mean()\n",
118+
" }\n",
119+
"\n",
120+
" # Add the results for the current cycle step size to the main dictionary\n",
121+
" all_results[\"results\"][i_cycle_step_size] = final_cycle_values\n",
122+
"\n",
123+
"# Save all results to a single JSON file in the tests/expected_outputs directory\n",
124+
"with open('../../tests/expected_outputs/KorakianitisMixedModelPP_expected_output.json', 'w') as f:\n",
125+
" json.dump(all_results, f, indent=4)"
126+
]
127+
}
128+
],
129+
"metadata": {
130+
"kernelspec": {
131+
"display_name": "venv",
132+
"language": "python",
133+
"name": "python3"
134+
},
135+
"language_info": {
136+
"codemirror_mode": {
137+
"name": "ipython",
138+
"version": 3
139+
},
140+
"file_extension": ".py",
141+
"mimetype": "text/x-python",
142+
"name": "python",
143+
"nbconvert_exporter": "python",
144+
"pygments_lexer": "ipython3",
145+
"version": "3.12.4"
146+
}
147+
},
148+
"nbformat": 4,
149+
"nbformat_minor": 2
150+
}

‎sandbox/notebooks/generate_KMM_KGO.ipynb

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -28,7 +28,7 @@
2828
},
2929
{
3030
"cell_type": "code",
31-
"execution_count": null,
31+
"execution_count": 2,
3232
"metadata": {},
3333
"outputs": [
3434
{
@@ -79,8 +79,6 @@
7979
" # Initializing the solver\n",
8080
" solver = Solver(model=model)\n",
8181
"\n",
82-
" print(f\"Running simulation with cycle step size: {i_cycle_step_size}\")\n",
83-
"\n",
8482
" # Initializing the parameter object\n",
8583
" parobj = KorakianitisMixedModel_parameters()\n",
8684
"\n",

‎src/ModularCirc/Components/HC_mixed_elastance_pp.py

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
from .ComponentBase import ComponentBase
2+
from ..HelperRoutines import activation_function_1, \
3+
chamber_volume_rate_change, \
4+
time_shift
5+
from ..Time import TimeClass
6+
7+
import pandas as pd
8+
import numpy as np
9+
10+
11+
def dvdt(t, q_in=None, q_out=None, v=None, v_ref=0.0, y=None):
12+
if y is not None:
13+
q_in, q_out, v = y[:3]
14+
if q_in - q_out > 0.0 or v > v_ref:
15+
return q_in - q_out
16+
else:
17+
return 0.0
18+
19+
def gen_active_p(E_act, v_ref):
20+
def active_p(v):
21+
return E_act * (v - v_ref)
22+
return active_p
23+
24+
def gen_active_dpdt(E_act):
25+
def active_dpdt(q_i, q_o):
26+
return E_act * (q_i - q_o)
27+
return active_dpdt
28+
29+
def gen_passive_p(E_pas, k_pas, v_ref):
30+
def passive_p(v):
31+
return E_pas * (np.exp(k_pas * (v - v_ref)) - 1.0)
32+
return passive_p
33+
34+
def gen_passive_dpdt(E_pas, k_pas, v_ref):
35+
def passive_dpdt(v, q_i, q_o):
36+
return E_pas * k_pas * np.exp(k_pas * (v - v_ref)) * (q_i - q_o)
37+
return passive_dpdt
38+
39+
def gen_total_p(_af, active_p, passive_p):
40+
def total_p(t, y):
41+
_af_t = _af(t)
42+
return _af_t * active_p(y) + passive_p(y)
43+
return total_p
44+
45+
def gen_total_dpdt(active_p, passive_p, _af, active_dpdt, passive_dpdt):
46+
def total_dpdt(t, y):
47+
_af_t = _af(t)
48+
_d_af_dt = _af(t, dt=True)
49+
return (_d_af_dt * active_p(y[0]) + _af_t * active_dpdt(y[1], y[2]) + passive_dpdt(y[0], y[1], y[2]))
50+
return total_dpdt
51+
52+
def gen_comp_v(E_pas, v_ref, k_pas):
53+
def comp_v(t, y):
54+
return v_ref + np.log(y[0] / E_pas + 1.0) / k_pas
55+
return comp_v
56+
57+
def gen_time_shifter(delay_, T):
58+
def time_shifter(t):
59+
return time_shift(t, delay_, T)
60+
return time_shifter
61+
62+
def gen__af(af, time_shifter, kwargs):
63+
varnames = [name for name in af.__code__.co_varnames if name != 'coeff' and name != 't']
64+
kwargs2 = {key: val for key,val in kwargs.items() if key in varnames}
65+
def _af(t, dt=False):
66+
return af(time_shifter(t), dt=dt, **kwargs2)
67+
return _af
68+
69+
class HC_mixed_elastance_pp(ComponentBase):
70+
def __init__(self,
71+
name:str,
72+
time_object: TimeClass,
73+
E_pas: float,
74+
E_act: float,
75+
k_pas: float,
76+
v_ref: float,
77+
v : float = None,
78+
p : float = None,
79+
af = activation_function_1,
80+
*args, **kwargs
81+
) -> None:
82+
super().__init__(name=name, time_object=time_object, v=v, p=p)
83+
self.E_pas = E_pas
84+
self.k_pas = k_pas
85+
self.E_act = E_act
86+
self.v_ref = v_ref
87+
self.eps = 1.0e-3
88+
self.kwargs = kwargs
89+
self.af = af
90+
91+
self.make_unique_io_state_variable(p_flag=True, q_flag=False)
92+
93+
@property
94+
def P(self):
95+
return self._P_i._u
96+
97+
def setup(self) -> None:
98+
E_pas = self.E_pas
99+
k_pas = self.k_pas
100+
E_act = self.E_act
101+
v_ref = self.v_ref
102+
eps = self.eps
103+
kwargs= self.kwargs
104+
T = self._to.tcycle
105+
af = self.af
106+
107+
time_shifter = gen_time_shifter(delay_=kwargs['delay'], T=T)
108+
_af = gen__af(af=af, time_shifter=time_shifter, kwargs=kwargs)
109+
110+
active_p = gen_active_p(E_act=E_act, v_ref=v_ref)
111+
active_dpdt = gen_active_dpdt(E_act=E_act)
112+
passive_p = gen_passive_p(E_pas=E_pas, k_pas=k_pas, v_ref=v_ref)
113+
passive_dpdt = gen_passive_dpdt(E_pas=E_pas, k_pas=k_pas, v_ref=v_ref)
114+
total_p = gen_total_p(_af=_af, active_p=active_p, passive_p=passive_p)
115+
total_dpdt = gen_total_dpdt(active_p=active_p, passive_p=passive_p,
116+
_af=_af, active_dpdt=active_dpdt, passive_dpdt=passive_dpdt)
117+
comp_v = gen_comp_v(E_pas=E_pas, v_ref=v_ref, k_pas=k_pas)
118+
119+
self._V.set_dudt_func(chamber_volume_rate_change,
120+
function_name='chamber_volume_rate_change')
121+
self._V.set_inputs(pd.Series({'q_in' :self._Q_i.name,
122+
'q_out':self._Q_o.name}))
123+
124+
self._P_i.set_dudt_func(total_dpdt, function_name='total_dpdt')
125+
self._P_i.set_inputs(pd.Series({'v' :self._V.name,
126+
'q_i':self._Q_i.name,
127+
'q_o':self._Q_o.name}))
128+
if self.p0 is None or self.p0 is np.NaN:
129+
self._P_i.set_i_func(total_p, function_name='total_p')
130+
self._P_i.set_i_inputs(pd.Series({'v':self._V.name}))
131+
else:
132+
self.P_i.loc[0] = self.p0
133+
if self.v0 is None or self.v0 is np.NaN:
134+
self._V.set_i_func(comp_v, function_name='comp_v')
135+
self._V.set_i_inputs(pd.Series({'p':self._P_i.name}))
136+
if (self.v0 is None or self.v0 is np.NaN) and (self.p0 is None or self.p0 is np.NaN):
137+
raise Exception("Solver needs at least the initial volume or pressure to be defined!")

‎src/ModularCirc/Components/__init__.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
from .ComponentBase import ComponentBase
22
from .HC_constant_elastance import HC_constant_elastance
33
from .HC_mixed_elastance import HC_mixed_elastance
4+
from .HC_mixed_elastance_pp import HC_mixed_elastance_pp
45
from .R_component import R_component
56
from .Rc_component import Rc_component
67
from .Rlc_component import Rlc_component

0 commit comments

Comments
 (0)