Adding a variation of KorakianitisMixedModel with passive component

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 (<V>), pressure (<P_i>) and flow (<Q_i>)\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
+}

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",

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!")

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