diff --git a/sandbox/notebooks/generate_KMMPP_KGO.ipynb b/sandbox/notebooks/generate_KMMPP_KGO.ipynb new file mode 100644 index 0000000..a088f16 --- /dev/null +++ b/sandbox/notebooks/generate_KMMPP_KGO.ipynb @@ -0,0 +1,150 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Description\n", + "\n", + "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", + "\n", + "With different step sizes, there are slightly different solutions reached, so tests need these different KGOs to have as comparisons." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import json\n", + "\n", + "from ModularCirc.Solver import Solver\n", + "\n", + "from ModularCirc.Models.KorakianitisMixedModelPP import KorakianitisMixedModelPP\n", + "from ModularCirc.Models.KorakianitisMixedModel_parameters import KorakianitisMixedModel_parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running simulation with cycle step size: 1\n", + "Running simulation with cycle step size: 1\n", + "Running simulation with cycle step size: 3\n", + "Running simulation with cycle step size: 3\n", + "Running simulation with cycle step size: 5\n", + "Running simulation with cycle step size: 5\n", + "Running simulation with cycle step size: 7\n", + "Running simulation with cycle step size: 7\n" + ] + } + ], + "source": [ + "time_setup_dict = {\n", + " 'name': 'TimeTest',\n", + " 'ncycles': 40,\n", + " 'tcycle': 1.0,\n", + " 'dt': 0.001,\n", + " 'export_min': 1\n", + " }\n", + "\n", + "# Start the loop with different cycle step size to generate JSON files for each loop:\n", + "cycle_step_sizes = [1, 3, 5, 7]\n", + "\n", + "# Dictionary to store the results for all cycle step sizes\n", + "all_results = {\n", + " \"metadata\": {\n", + " \"description\": \"Results for different cycle_step_size values\",\n", + " \"cycle_step_sizes\": cycle_step_sizes\n", + " },\n", + " \"results\": {}\n", + "}\n", + "\n", + "# Loop through each cycle step size\n", + "for i_cycle_step_size in cycle_step_sizes:\n", + "\n", + " print(f\"Running simulation with cycle step size: {i_cycle_step_size}\")\n", + "\n", + " # Initializing the parameter object\n", + " parobj = KorakianitisMixedModel_parameters()\n", + "\n", + " # Initializing the model \n", + " model = KorakianitisMixedModelPP(time_setup_dict=time_setup_dict, \n", + " parobj=parobj, \n", + " suppress_printing=True)\n", + "\n", + " # Initializing the solver\n", + " solver = Solver(model=model)\n", + "\n", + " # Initializing the parameter object\n", + " parobj = KorakianitisMixedModel_parameters()\n", + "\n", + " # Initializing the model \n", + " model = KorakianitisMixedModelPP(time_setup_dict=time_setup_dict, \n", + " parobj=parobj, \n", + " suppress_printing=True)\n", + "\n", + " # Initializing the solver\n", + " solver = Solver(model=model)\n", + "\n", + " # Solver is being setup: switching off console printing and setting the solver method to \"LSODA\"\n", + " solver.setup(suppress_output=True, \n", + " method='LSODA',\n", + " step=i_cycle_step_size)\n", + "\n", + " # Running the model\n", + " solver.solve()\n", + "\n", + " # Define the indexes of the equivalent to the last cycles\n", + " tind_fin = np.arange(start=model.time_object.n_t-model.time_object.n_c,\n", + " stop=model.time_object.n_t)\n", + "\n", + " # Dictionary to store the final cycle values for each component for the current step size\n", + " final_cycle_values = {}\n", + "\n", + " # From each of the components, retrieve the volume (), pressure () and flow ()\n", + " for key, value in model.components.items():\n", + " final_cycle_values[key] = {\n", + " 'V': value.V.values[tind_fin].mean(),\n", + " 'P_i': value.P_i.values[tind_fin].mean(),\n", + " 'Q_i': value.Q_i.values[tind_fin].mean()\n", + " }\n", + "\n", + " # Add the results for the current cycle step size to the main dictionary\n", + " all_results[\"results\"][i_cycle_step_size] = final_cycle_values\n", + "\n", + "# Save all results to a single JSON file in the tests/expected_outputs directory\n", + "with open('../../tests/expected_outputs/KorakianitisMixedModelPP_expected_output.json', 'w') as f:\n", + " json.dump(all_results, f, indent=4)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/sandbox/notebooks/generate_KMM_KGO.ipynb b/sandbox/notebooks/generate_KMM_KGO.ipynb index 3fe5986..1aa9b5e 100644 --- a/sandbox/notebooks/generate_KMM_KGO.ipynb +++ b/sandbox/notebooks/generate_KMM_KGO.ipynb @@ -28,7 +28,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -79,8 +79,6 @@ " # Initializing the solver\n", " solver = Solver(model=model)\n", "\n", - " print(f\"Running simulation with cycle step size: {i_cycle_step_size}\")\n", - "\n", " # Initializing the parameter object\n", " parobj = KorakianitisMixedModel_parameters()\n", "\n", diff --git a/src/ModularCirc/Components/HC_mixed_elastance_pp.py b/src/ModularCirc/Components/HC_mixed_elastance_pp.py new file mode 100644 index 0000000..f90680e --- /dev/null +++ b/src/ModularCirc/Components/HC_mixed_elastance_pp.py @@ -0,0 +1,137 @@ +from .ComponentBase import ComponentBase +from ..HelperRoutines import activation_function_1, \ + chamber_volume_rate_change, \ + time_shift +from ..Time import TimeClass + +import pandas as pd +import numpy as np + + +def dvdt(t, q_in=None, q_out=None, v=None, v_ref=0.0, y=None): + if y is not None: + q_in, q_out, v = y[:3] + if q_in - q_out > 0.0 or v > v_ref: + return q_in - q_out + else: + return 0.0 + +def gen_active_p(E_act, v_ref): + def active_p(v): + return E_act * (v - v_ref) + return active_p + +def gen_active_dpdt(E_act): + def active_dpdt(q_i, q_o): + return E_act * (q_i - q_o) + return active_dpdt + +def gen_passive_p(E_pas, k_pas, v_ref): + def passive_p(v): + return E_pas * (np.exp(k_pas * (v - v_ref)) - 1.0) + return passive_p + +def gen_passive_dpdt(E_pas, k_pas, v_ref): + def passive_dpdt(v, q_i, q_o): + return E_pas * k_pas * np.exp(k_pas * (v - v_ref)) * (q_i - q_o) + return passive_dpdt + +def gen_total_p(_af, active_p, passive_p): + def total_p(t, y): + _af_t = _af(t) + return _af_t * active_p(y) + passive_p(y) + return total_p + +def gen_total_dpdt(active_p, passive_p, _af, active_dpdt, passive_dpdt): + def total_dpdt(t, y): + _af_t = _af(t) + _d_af_dt = _af(t, dt=True) + return (_d_af_dt * active_p(y[0]) + _af_t * active_dpdt(y[1], y[2]) + passive_dpdt(y[0], y[1], y[2])) + return total_dpdt + +def gen_comp_v(E_pas, v_ref, k_pas): + def comp_v(t, y): + return v_ref + np.log(y[0] / E_pas + 1.0) / k_pas + return comp_v + +def gen_time_shifter(delay_, T): + def time_shifter(t): + return time_shift(t, delay_, T) + return time_shifter + +def gen__af(af, time_shifter, kwargs): + varnames = [name for name in af.__code__.co_varnames if name != 'coeff' and name != 't'] + kwargs2 = {key: val for key,val in kwargs.items() if key in varnames} + def _af(t, dt=False): + return af(time_shifter(t), dt=dt, **kwargs2) + return _af + +class HC_mixed_elastance_pp(ComponentBase): + def __init__(self, + name:str, + time_object: TimeClass, + E_pas: float, + E_act: float, + k_pas: float, + v_ref: float, + v : float = None, + p : float = None, + af = activation_function_1, + *args, **kwargs + ) -> None: + super().__init__(name=name, time_object=time_object, v=v, p=p) + self.E_pas = E_pas + self.k_pas = k_pas + self.E_act = E_act + self.v_ref = v_ref + self.eps = 1.0e-3 + self.kwargs = kwargs + self.af = af + + self.make_unique_io_state_variable(p_flag=True, q_flag=False) + + @property + def P(self): + return self._P_i._u + + def setup(self) -> None: + E_pas = self.E_pas + k_pas = self.k_pas + E_act = self.E_act + v_ref = self.v_ref + eps = self.eps + kwargs= self.kwargs + T = self._to.tcycle + af = self.af + + time_shifter = gen_time_shifter(delay_=kwargs['delay'], T=T) + _af = gen__af(af=af, time_shifter=time_shifter, kwargs=kwargs) + + active_p = gen_active_p(E_act=E_act, v_ref=v_ref) + active_dpdt = gen_active_dpdt(E_act=E_act) + passive_p = gen_passive_p(E_pas=E_pas, k_pas=k_pas, v_ref=v_ref) + passive_dpdt = gen_passive_dpdt(E_pas=E_pas, k_pas=k_pas, v_ref=v_ref) + total_p = gen_total_p(_af=_af, active_p=active_p, passive_p=passive_p) + total_dpdt = gen_total_dpdt(active_p=active_p, passive_p=passive_p, + _af=_af, active_dpdt=active_dpdt, passive_dpdt=passive_dpdt) + comp_v = gen_comp_v(E_pas=E_pas, v_ref=v_ref, k_pas=k_pas) + + self._V.set_dudt_func(chamber_volume_rate_change, + function_name='chamber_volume_rate_change') + self._V.set_inputs(pd.Series({'q_in' :self._Q_i.name, + 'q_out':self._Q_o.name})) + + self._P_i.set_dudt_func(total_dpdt, function_name='total_dpdt') + self._P_i.set_inputs(pd.Series({'v' :self._V.name, + 'q_i':self._Q_i.name, + 'q_o':self._Q_o.name})) + if self.p0 is None or self.p0 is np.NaN: + self._P_i.set_i_func(total_p, function_name='total_p') + self._P_i.set_i_inputs(pd.Series({'v':self._V.name})) + else: + self.P_i.loc[0] = self.p0 + if self.v0 is None or self.v0 is np.NaN: + self._V.set_i_func(comp_v, function_name='comp_v') + self._V.set_i_inputs(pd.Series({'p':self._P_i.name})) + if (self.v0 is None or self.v0 is np.NaN) and (self.p0 is None or self.p0 is np.NaN): + raise Exception("Solver needs at least the initial volume or pressure to be defined!") diff --git a/src/ModularCirc/Components/__init__.py b/src/ModularCirc/Components/__init__.py index c79b8d9..ff9fa56 100644 --- a/src/ModularCirc/Components/__init__.py +++ b/src/ModularCirc/Components/__init__.py @@ -1,6 +1,7 @@ from .ComponentBase import ComponentBase from .HC_constant_elastance import HC_constant_elastance from .HC_mixed_elastance import HC_mixed_elastance +from .HC_mixed_elastance_pp import HC_mixed_elastance_pp from .R_component import R_component from .Rc_component import Rc_component from .Rlc_component import Rlc_component diff --git a/src/ModularCirc/Models/KorakianitisMixedModelPP.py b/src/ModularCirc/Models/KorakianitisMixedModelPP.py new file mode 100644 index 0000000..fbf946c --- /dev/null +++ b/src/ModularCirc/Models/KorakianitisMixedModelPP.py @@ -0,0 +1,114 @@ +from .OdeModel import OdeModel +from .KorakianitisMixedModel_parameters import KorakianitisMixedModel_parameters, TIME_SETUP_DICT +from .KorakianitisMixedModel_parameters import TIME_SETUP_DICT as TEMPLATE_TIME_SETUP_DICT +from .ParametersObject import ParametersObject as po +from ..Components import Rlc_component, Valve_simple_bernoulli, HC_mixed_elastance_pp + +FULL_NAMES =[ + 'LeftA', + 'MiValve', + 'LeftV', + 'AoV', + 'SysAoSin', + 'SysArt', + 'SysVen', + 'RightA', + 'TriValve', + 'RightV', + 'PulV', + 'PulArtSin', + 'PulArt', + 'PulVen', +] + +class KorakianitisMixedModelPP(OdeModel): + def __init__(self, time_setup_dict, parobj:po=KorakianitisMixedModel_parameters, suppress_printing:bool=False) -> None: + super().__init__(time_setup_dict) + self.name = 'KorakianitisModel' + + if not suppress_printing: print(parobj) + + # The components... + for key, name in zip(parobj.components.keys(), FULL_NAMES): + if key in parobj._vessels: + class_ = Rlc_component + elif key in parobj._valves: + class_ = Valve_simple_bernoulli + elif key in parobj._chambers: + class_ = HC_mixed_elastance_pp + else: + raise Exception(f'Component name {key} not in the model list.') + self.components[key] = class_(name=name, + time_object=self.time_object, + **parobj[key].to_dict()) + if key not in parobj._valves: + self.set_v_sv(key) + # else: + # self.set_phi_sv(key) + self.components[key].setup() + + self.connect_modules(self.components['lv'], + self.components['ao'], + plabel='p_lv', + qlabel='q_ao') + self.connect_modules(self.components['ao'], + self.components['sas'], + plabel='p_sas', + qlabel='q_ao') + self.connect_modules(self.components['sas'], + self.components['sat'], + plabel='p_sat', + qlabel='q_sas') + self.connect_modules(self.components['sat'], + self.components['svn'], + plabel='p_svn', + qlabel='q_sat') + self.connect_modules(self.components['svn'], + self.components['ra'], + plabel='p_ra', + qlabel='q_svn') + self.connect_modules(self.components['ra'], + self.components['ti'], + plabel='p_ra', + qlabel='q_ti') + self.connect_modules(self.components['ti'], + self.components['rv'], + plabel='p_rv', + qlabel='q_ti') + self.connect_modules(self.components['rv'], + self.components['po'], + plabel='p_rv', + qlabel='q_po') + self.connect_modules(self.components['po'], + self.components['pas'], + plabel='p_pas', + qlabel='q_po') + self.connect_modules(self.components['pas'], + self.components['pat'], + plabel='p_pat', + qlabel='q_pas') + self.connect_modules(self.components['pat'], + self.components['pvn'], + plabel='p_pvn', + qlabel='q_pat') + self.connect_modules(self.components['pvn'], + self.components['la'], + plabel='p_la', + qlabel='q_pvn') + self.connect_modules(self.components['la'], + self.components['mi'], + plabel='p_la', + qlabel='q_mi') + self.connect_modules(self.components['mi'], + self.components['lv'], + plabel='p_lv', + qlabel='q_mi') + + for component in self.components.values(): + component.setup() + + # def set_phi_sv(self, comp_key:str) -> None: + # phi_key = 'phi_' + comp_key + # self._state_variable_dict[phi_key] = self.components[comp_key]._PHI + # self._state_variable_dict[phi_key].set_name(phi_key) + # self.all_sv_data[phi_key] = self.components[comp_key].PHI diff --git a/tests/expected_outputs/KorakianitisMixedModelPP_expected_output.json b/tests/expected_outputs/KorakianitisMixedModelPP_expected_output.json new file mode 100644 index 0000000..fa2c22e --- /dev/null +++ b/tests/expected_outputs/KorakianitisMixedModelPP_expected_output.json @@ -0,0 +1,301 @@ +{ + "metadata": { + "description": "Results for different cycle_step_size values", + "cycle_step_sizes": [ + 1, + 3, + 5, + 7 + ] + }, + "results": { + "1": { + "la": { + "V": 125.43780896363097, + "P_i": 1.8918443840149175, + "Q_i": 85.63841933794393 + }, + "mi": { + "V": 0.0, + "P_i": 1.8918443840149175, + "Q_i": 85.12217263877217 + }, + "lv": { + "V": 78.43302167656424, + "P_i": 35.726596708394545, + "Q_i": 85.12217263877217 + }, + "ao": { + "V": 0.0, + "P_i": 35.726596708394545, + "Q_i": 85.46707071257764 + }, + "sas": { + "V": 7.935169985533782, + "P_i": 99.18962481915754, + "Q_i": 85.46707071257764 + }, + "sat": { + "V": 158.29328734571033, + "P_i": 98.9333045910714, + "Q_i": 85.46670858628629 + }, + "svn": { + "V": 146.52144876746553, + "P_i": 7.14738774475429, + "Q_i": 85.78154959253648 + }, + "ra": { + "V": 53.80419144699683, + "P_i": 0.7409492646472583, + "Q_i": 85.41917973476043 + }, + "ti": { + "V": 0.0, + "P_i": 0.7409492646472583, + "Q_i": 85.21441802598888 + }, + "rv": { + "V": 66.18665308761932, + "P_i": 11.556955202253105, + "Q_i": 85.21441802598888 + }, + "po": { + "V": 0.0, + "P_i": 11.556955202253105, + "Q_i": 84.76151384790603 + }, + "pas": { + "V": 5.188167663743965, + "P_i": 28.82315368748437, + "Q_i": 84.76151384790603 + }, + "pat": { + "V": 108.88391561186756, + "P_i": 28.653662003122296, + "Q_i": 84.75736435946384 + }, + "pvn": { + "V": 49.31633545087104, + "P_i": 2.4056749000425808, + "Q_i": 84.6698041042773 + } + }, + "3": { + "la": { + "V": 125.43671555581082, + "P_i": 1.8918743966439486, + "Q_i": 85.63899452725995 + }, + "mi": { + "V": 0.0, + "P_i": 1.8918743966439486, + "Q_i": 85.12276623473828 + }, + "lv": { + "V": 78.43370556654797, + "P_i": 35.726732080652454, + "Q_i": 85.12276623473828 + }, + "ao": { + "V": 0.0, + "P_i": 35.726732080652454, + "Q_i": 85.46775977332928 + }, + "sas": { + "V": 7.935200306088251, + "P_i": 99.19000382611331, + "Q_i": 85.46775977332928 + }, + "sat": { + "V": 158.29389016355236, + "P_i": 98.93368135221976, + "Q_i": 85.4672576752595 + }, + "svn": { + "V": 146.52117445767013, + "P_i": 7.14737436378872, + "Q_i": 85.7819140635977 + }, + "ra": { + "V": 53.802429334604184, + "P_i": 0.7409260081657724, + "Q_i": 85.41931140830596 + }, + "ti": { + "V": 0.0, + "P_i": 0.7409260081657724, + "Q_i": 85.21423254356965 + }, + "rv": { + "V": 66.18739846234556, + "P_i": 11.556977257422208, + "Q_i": 85.21423254356965 + }, + "po": { + "V": 0.0, + "P_i": 11.556977257422208, + "Q_i": 84.76185074673485 + }, + "pas": { + "V": 5.1881849730035094, + "P_i": 28.823249850020396, + "Q_i": 84.76185074673485 + }, + "pat": { + "V": 108.88427972232917, + "P_i": 28.653757821665934, + "Q_i": 84.75771185332971 + }, + "pvn": { + "V": 49.317021458053645, + "P_i": 2.4057083638075087, + "Q_i": 84.6700042863683 + } + }, + "5": { + "la": { + "V": 125.13647946978168, + "P_i": 1.885908356875621, + "Q_i": 85.85236009026312 + }, + "mi": { + "V": 0.0, + "P_i": 1.885908356875621, + "Q_i": 85.19668336916295 + }, + "lv": { + "V": 78.1973923686021, + "P_i": 35.56087846444398, + "Q_i": 85.19668336916295 + }, + "ao": { + "V": 0.0, + "P_i": 35.56087846444398, + "Q_i": 85.18418104976712 + }, + "sas": { + "V": 7.8979264354362435, + "P_i": 98.72408044296233, + "Q_i": 85.18418104976712 + }, + "sat": { + "V": 157.54983238781537, + "P_i": 98.46864524238467, + "Q_i": 85.16341090231866 + }, + "svn": { + "V": 146.6416237510697, + "P_i": 7.1532499390763835, + "Q_i": 85.3416490128491 + }, + "ra": { + "V": 54.05281190813093, + "P_i": 0.7448176773630025, + "Q_i": 85.44576348951175 + }, + "ti": { + "V": 0.0, + "P_i": 0.7448176773630025, + "Q_i": 85.36639423785674 + }, + "rv": { + "V": 66.49674115758373, + "P_i": 11.635142251196598, + "Q_i": 85.36639423785674 + }, + "po": { + "V": 0.0, + "P_i": 11.635142251196598, + "Q_i": 85.34931733799588 + }, + "pas": { + "V": 5.22157011712042, + "P_i": 29.008722872892957, + "Q_i": 85.34931733799588 + }, + "pat": { + "V": 109.58466079740198, + "P_i": 28.838068630896718, + "Q_i": 85.33840482829541 + }, + "pvn": { + "V": 49.220961607049915, + "P_i": 2.4010225174172, + "Q_i": 85.2801603466206 + } + }, + "7": { + "la": { + "V": 125.12976179137765, + "P_i": 1.8855774322008392, + "Q_i": 85.77008943682864 + }, + "mi": { + "V": 0.0, + "P_i": 1.8855774322008392, + "Q_i": 85.13276321844744 + }, + "lv": { + "V": 78.21261330483122, + "P_i": 35.5882402798702, + "Q_i": 85.13276321844744 + }, + "ao": { + "V": 0.0, + "P_i": 35.5882402798702, + "Q_i": 85.1574677529472 + }, + "sas": { + "V": 7.905076155825749, + "P_i": 98.81345194783403, + "Q_i": 85.1574677529472 + }, + "sat": { + "V": 157.69291523694326, + "P_i": 98.55807202308897, + "Q_i": 85.15256225938438 + }, + "svn": { + "V": 146.71642465198522, + "P_i": 7.156898763511504, + "Q_i": 85.42190218338914 + }, + "ra": { + "V": 54.040798347985394, + "P_i": 0.7446781355125455, + "Q_i": 85.49627503998612 + }, + "ti": { + "V": 0.0, + "P_i": 0.7446781355125455, + "Q_i": 85.3945765573658 + }, + "rv": { + "V": 66.45364929420494, + "P_i": 11.619722813932832, + "Q_i": 85.3945765573658 + }, + "po": { + "V": 0.0, + "P_i": 11.619722813932832, + "Q_i": 85.30189730432218 + }, + "pas": { + "V": 5.21424536571159, + "P_i": 28.968029809510313, + "Q_i": 85.30189730432218 + }, + "pat": { + "V": 109.4304574902965, + "P_i": 28.79748881323648, + "Q_i": 85.27987458268527 + }, + "pvn": { + "V": 49.20405836084357, + "P_i": 2.400197968821811, + "Q_i": 85.15154807247924 + } + } + } +} diff --git a/tests/test_KorakianitisMixedModelPP.py b/tests/test_KorakianitisMixedModelPP.py new file mode 100644 index 0000000..d472884 --- /dev/null +++ b/tests/test_KorakianitisMixedModelPP.py @@ -0,0 +1,171 @@ +import unittest +import numpy as np +import json +import os +import logging +from ModularCirc.Models.KorakianitisMixedModelPP import KorakianitisMixedModelPP +from ModularCirc.Models.KorakianitisMixedModel_parameters import KorakianitisMixedModel_parameters +from ModularCirc.Solver import Solver + +# Configure logging +logging.basicConfig(level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s') + +# Define global constants for tolerances +RELATIVE_TOLERANCE = 1e-3 + +class TestKorakianitisMixedModelPP(unittest.TestCase): + + def setUp(self): + """ + Set up the test environment for the KorakianitisMixedModelPP. + + This method initializes the necessary components and configurations + required for testing the KorakianitisMixedModelPP. It performs the following: + + - Sets a random seed for reproducibility. + - Defines the base directory for file paths. + - Configures the simulation time setup parameters, including the number + of cycles, cycle duration, time step size, and minimum export cycles. + - Initializes the parameter object for the KorakianitisMixedModelPP. + - Creates an instance of the KorakianitisMixedModelPP with the specified + time setup and parameter object, suppressing console output. + - Initializes the solver for the model and configures it with the + "LSODA" method, suppressing console output. + - Loads expected output values from a JSON file for validation during tests. + - Verifies the existence of the expected output file and raises an + assertion error if the file is not found. + + Raises: + AssertionError: If the expected output file is not found. + """ + + # Set a random seed for reproducibility + np.random.seed(42) + + # Define the base directory for file paths + self.base_dir = os.path.dirname(__file__) + + # Define the duration of the simulation (no of cycles), duration of the cycle, maximum time step size, and minimum number of cycles to run + self.time_setup_dict = { + 'name': 'TimeTest', + 'ncycles': 40, + 'tcycle': 1.0, + 'dt': 0.001, + 'export_min': 1 + } + + # Initializing the parameter object + self.parobj = KorakianitisMixedModel_parameters() + + # Initializing the model + self.model = KorakianitisMixedModelPP(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, + suppress_printing=True) + + # Initializing the solver + self.solver = Solver(model=self.model) + + # Solver is being setup: switching off console printing and setting the solver method to "LSODA" + self.solver.setup(suppress_output=True, method='LSODA',step=1) + + # Load expected values from a JSON file + output_file_path = os.path.join(self.base_dir, 'expected_outputs', 'KorakianitisMixedModelPP_expected_output.json') + + # Verify the file exists + self.assertTrue(os.path.exists(output_file_path), f"Expected output file not found: {output_file_path}") + + with open(output_file_path, 'r') as f: + self.expected_values = json.load(f) + + def test_model_initialization(self): + ''' + Testing the initialization of the solver, model and parameter objects. + ''' + # Verify model is an instance of + self.assertIsInstance(self.model, KorakianitisMixedModelPP) + # Verify model has attribute + self.assertTrue(hasattr(self.solver.model, 'components')) + # Verify is a component + self.assertIn('lv', self.solver.model.components) + # Verify correct assignment of parameters from parobj to model + self.assertEqual(self.solver.model.components['lv'].E_pas, self.parobj.components['lv']['E_pas']) + self.assertEqual(self.solver.model.components['ao'].CQ, self.parobj.components['ao']['CQ']) + + def test_solver_initialization(self): + # Verify is an instance of + self.assertIsInstance(self.solver, Solver) + # Verify the instance of that is an atribute of is the same as the original + self.assertEqual(self.solver.model, self.model) + + def test_solver_run(self): + """ + Test the functionality of the solver by running it with different step sizes and verifying its behavior. + This test ensures that: + 1. The solver can be configured and run with various step sizes. + 2. The solver converges successfully for each step size. + 3. The state variables of the model components are updated correctly after solving. + 4. The computed results match the expected values within a specified tolerance. + + """ + + cycle_step_sizes = [1, 3, 5, 7] # Define the step sizes to test + + for i_cycle_step_size in cycle_step_sizes: + + # Use logging to print the current step size + logging.info(f"Testing solver with step size: {i_cycle_step_size}") + + with self.subTest(cycle_step_size=i_cycle_step_size): + + # Initializing the parameter object + self.parobj = KorakianitisMixedModel_parameters() + + # Initializing the model + self.model = KorakianitisMixedModelPP(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, + suppress_printing=True) + + # Initializing the solver + self.solver = Solver(model=self.model) + + # Reconfigure the solver with the current step size + self.solver.setup(suppress_output=True, method='LSODA', step=i_cycle_step_size) + + # Running the model + self.solver.solve() + + # Verify the solver converged + self.assertTrue(self.solver.converged or self.solver._Nconv is not None) + + # Verifying the model changed the state variables stored within components. + self.assertTrue(len(self.solver.model.components['lv'].V.values) > 0) + self.assertTrue(len(self.solver.model.components['lv'].P_i.values) > 0) + + # Redefine tind based on how many heart cycle have actually been necessary to reach steady state + self.tind_fin = np.arange(start=self.model.time_object.n_t-self.model.time_object.n_c, + stop=(self.model.time_object.n_t)) + + # Retrieve the component state variables, compute the mean of the values during the last cycle and store them within + # the new solution dictionary + new_dict = {} + for key, value in self.model.components.items(): + + new_dict[key] = { + 'V': value.V.values[self.tind_fin].mean(), + 'P_i': value.P_i.values[self.tind_fin].mean(), + 'Q_i': value.Q_i.values[self.tind_fin].mean() + } + + # Check that the values are the same as the expected values + expected_ndarray = np.array( + [self.expected_values["results"][str(i_cycle_step_size)][key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()] + ) + new_ndarray = np.array([new_dict[key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()]) + test_ndarray = np.where(np.abs(expected_ndarray) > 1e-6, + np.abs((expected_ndarray - new_ndarray) / expected_ndarray), + np.abs((expected_ndarray - new_ndarray))) + self.assertTrue((test_ndarray < RELATIVE_TOLERANCE).all()) + + +if __name__ == '__main__': + unittest.main()