Build your own model class

Todo

Explain how to create a custom model class

When the model is created, it must be parameterized properly. It is explained here how to create your own reference model state.

The model VanOsta2023 is build following this standard. This model is constructed by the following code.

  1"""
  2VanOsta2023 model.
  3
  4"""
  5
  6from circadapt.model import Model
  7from circadapt.adapt import ModelAdapt
  8from circadapt.plot import triseg2022
  9import matplotlib.pyplot as plt
 10import numpy as np
 11
 12
 13class VanOsta2023(Model, ModelAdapt):
 14    def __init__(self,
 15                 solver: str = None,
 16                 path_to_circadapt: str = None,
 17                 model_state: dict = None,
 18                 ):
 19        if solver is None:
 20            solver = 'backward_differential'
 21
 22        self._local_save_reference = True
 23
 24        ModelAdapt.__init__(self)
 25        Model.__init__(self,
 26                       solver,
 27                       path_to_circadapt=path_to_circadapt,
 28                       model_state=model_state,
 29                       )
 30
 31    def _update_settings(self):
 32        # rename model components for easy input/output
 33        self._module_rename['Wall2022'] = 'Wall'
 34        self._module_rename['Patch2022'] = 'Patch'
 35        self._module_rename['Chamber2022'] = 'Chamber'
 36        self._module_rename['Valve2022'] = 'Valve'
 37        self._module_rename['TriSeg2022'] = 'TriSeg'
 38        self._module_rename['PressureFlowControl'] = 'PFC'
 39
 40    def build(self):
 41        # Circulation
 42        self.add_smart_component('ArtVen', build='SystemicCirculation')
 43        self.add_smart_component('ArtVen', build='PulmonaryCirculation')
 44        self.add_smart_component('Heart', patch_type='Patch2022')
 45        self.add_smart_component('Timings')
 46        self.add_smart_component('PressureFlowControl')
 47
 48        # manually set papillary muscles
 49        self.set_component("Peri.RaRv.wPapMus", "Peri.TriSeg.wRv")
 50        self.set_component("Peri.LaLv.wPapMus", "Peri.TriSeg.wLv")
 51
 52    def set_reference(self):
 53        self['Chamber']['buckling'] = False
 54        self['Valve']['soft_closure'] = True
 55        self['Valve']['papillary_muscles'] = True
 56
 57        self['Solver']['dt'] = 0.001
 58        self['Solver']['dt_export'] = 0.002
 59        self.set('Solver.order',  2)
 60
 61        self.set('Model.t_cycle', 0.85)
 62
 63        self.set('Model.PFC.fac',  0.5)
 64        self.set('Model.PFC.epsilon', 0.1)
 65        self.set('Model.PFC.fac_pfc', 1)
 66        self.set('Model.PFC.stable_threshold', 0.001)
 67
 68        # Set parameters Circulation
 69        self['Tube0D']['A0'] = [0.00049828, 0.00049899, 0.00046929, 0.00050574]
 70        self['Tube0D']['A_wall'] = [1.13187792e-04, 3.72289316e-05, 7.27489286e-05, 4.31617949e-05]
 71        self['Tube0D']['k'] = [1.66666667, 2.33333333, 1.66666667, 2.33333333]
 72        self['Tube0D']['l'] = [0.4, 0.4, 0.19, 0.19]
 73        self['Tube0D']['p0'] = [12155.25602777,   297.44066943,  1912.6684696 ,   580.51425824]
 74
 75        self['ArtVen']['p0'] = [6299.86481149, 1000.        ]
 76        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
 77        self['ArtVen']['k'] = [1.  , 1.88]
 78
 79        # Set volume state variables
 80        # not needed to reset, only if initial model crashes
 81        if True:
 82            self.set('Model.SyArt.V', 2e-4)
 83            self.set('Model.PuArt.V', 1e-4)
 84            self.set('Model.SyVen.V', 4e-4)
 85            self.set('Model.PuVen.V', 2e-4)
 86
 87            self.set('Model.Peri.TriSeg.V', 56e-6)
 88            self.set('Model.Peri.TriSeg.Y', 36e-3)
 89
 90            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
 91            self.set('Model.Peri.TriSeg.cRv.V',  110e-6)
 92            self.set('Model.Peri.La.V',  25e-6)
 93            self.set('Model.Peri.Ra.V',  150e-6)
 94
 95        # all wall parameters
 96        self['Patch']['l_se0'] = 0.04
 97        self['Patch']['l_s_ref'] = 2.0
 98        self['Patch']['l_s0'] = [1.79, 1.79, 1.8 , 1.8 , 1.8 ]
 99        self['Patch']['dl_s_pas'] = 0.6
100        self['Patch']['k1'] = 10
101        self['Patch']['dt'] = 0
102        self['Patch']['C_rest'] = 0
103        self['Patch']['l_si0'] = 1.51
104        self['Patch']['LDAD'] = [1.057, 1.057, 0.74 , 0.74 , 0.74 ]
105        self['Patch']['ADO'] = [0.081, 0.081, 0.75 , 0.75 , 0.75 ]
106        self['Patch']['LDCC'] = [4. , 4. , 3., 3., 3.]
107        self['Patch']['v_max'] = 7.
108        self['Patch']['v_max'][:2] = 14.
109
110        self['Patch']['tr'] = [0.4 , 0.4 , 0.21, 0.21, 0.21]
111        self['Patch']['td'] = [0.4 , 0.4 , 0.24, 0.24, 0.24]
112
113        # wall specific
114        self['Patch']['Sf_act'] = [100000., 100000., 120000., 120000., 120000.]
115        self['Patch']['Am_ref'] = [0.00400721, 0.00311263, 0.00918281, 0.00518683, 0.01321999]
116        self['Patch']['V_wall'] = [1.59935455e-05, 7.72140495e-06, 1.16441768e-04, 4.51241560e-05,
117               6.88196851e-05]
118        self['Patch']['Sf_pas'] = [  9.51064014,  10.3647024 , 424.3537027 , 417.74300178,
119               439.21574206]
120
121        # Valves
122        self['Valve']['adaptation_A_open_fac'] = [1.  , 1.12, 1.  , 1.  , 1.12, 1.  ]
123        self['Valve']['A_open'] = [0.00050044, 0.00052598, 0.00046962, 0.00050745, 0.0005583 ,
124               0.00049849]
125        self['Valve']['A_leak'] = [2.64705882e-04, 2.64705882e-10, 2.64705882e-10, 2.64705882e-04,
126               2.64705882e-10, 2.64705882e-10]
127        self['Valve']['l'] = 0.01626978
128        self['Valve']['rho_b'] = 1050
129        self['Valve']['papillary_muscles'] = True
130        self['Valve']['papillary_muscles_slope'] = 100
131        self['Valve']['papillary_muscles_min'] = 0.1
132        self['Valve']['papillary_muscles_A_open_fac'] = 0.1
133        self['Valve']['soft_closure'] = True
134
135        self['Bag']['k'] = 10
136        self['Bag']['V_ref'] = [0.00057937]
137        self['Bag']['p_ref'] = 100
138
139        self.set('Model.Peri.TriSeg.Y', 0.035)
140
141        self['Timings']['law_tau_av'] = 2
142        self['Timings']['c_tau_av0'] = 0.189
143        self['Timings']['c_tau_av1'] = -0.029
144
145        # Adaptation
146        self['Patch']['adapt_gamma'] = 0.5
147        self.set('Solver.store_beats', 1)
148
149        self['PFC']['stable_threshold'] = 1e-3
150
151        self['Patch']['SfPasMaxT'] = [3600., 3600.,  4200.,  4200.,  4200.]
152        self['Patch']['FacSfActT'] = [0.28, 0.28, 0.69, 0.69, 0.69]
153        self['Patch']['SfPasActT'] = [2800., 2800., 6600., 6600., 6600.]
154        self['Patch']['LsPasActT'] = [3.  , 3.  , 2.31, 2.31, 2.31]
155
156        self.run(stable=True)
157
158        return
159
160
161    def get_unittest_targets(self):
162        """Hardcoded results after initializing and running 1 beat."""
163        return {
164            'LVEDV': 130.9,
165            'LVESV':  58.5,
166            }
167
168    def get_unittest_results(self, model):
169        """Real-time results after initializing and running 1 beat."""
170        LVEDV = np.max(model['Cavity']['V'][:, 'cLv'])*1e6
171        LVESV = np.min(model['Cavity']['V'][:, 'cLv'])*1e6
172        return {
173            'LVEDV': LVEDV,
174            'LVESV': LVESV,
175            }
176
177    def plot(self, fig=None):
178        # TODO
179        self.plot_extended(fig)
180
181    def plot_extended(self, fig=None):
182        if fig is None:
183            fig = 1
184        if isinstance(fig, int):
185            fig = plt.figure(fig, clear=True, figsize=(12, 8))
186
187        # Settings
188        grid_size = [32, 32]
189
190        def get_lim(module, signal, locs=slice(None, None, None)):
191            signal = self[module][signal][:, locs]
192            lim = np.array([np.min(signal), np.max(signal)])
193            lim += np.array([-1, 1]) * 0.1*np.diff(lim)
194            return lim
195
196        lim_V = get_lim('Cavity', 'V', ['cRv', 'Ra', 'La', 'cLv']) * 1e6
197        lim_V[0] = 0
198        lim_p = get_lim('Cavity', 'p', ['cRv', 'Ra', 'La', 'cLv']) / 133
199        lim_p[0] = np.min([lim_p[0], 0])
200        lim_Ls = get_lim('Patch', 'l_s')
201        lim_Sf = get_lim('Patch', 'Sf') * 1e-3
202        lim_q = get_lim('Valve', 'q',
203                        ['LaLv', 'RaRv', 'LvSyArt', 'RaPuArt']) * 1e6
204
205        all_lim = [lim_V, lim_p, lim_Ls, lim_Sf, lim_q]
206        if (np.any(np.isnan(all_lim)) or np.any(np.isinf(all_lim))):
207            lim_V = [0, 200]
208            lim_p = [0, 150]
209            lim_Ls = [1.5, 2.0]
210            lim_Sf = [0, 100]
211            lim_q = [-1e-3, 1e-3]
212
213        # Pressure Volume plot
214        axPV = plt.subplot2grid(grid_size, (0, 17), rowspan=15, colspan=15, fig=fig)
215        axPV.plot(self['Cavity']['V'][:, 'cLv']*1e6, self['Cavity']['p'][:, 'cLv']/133)
216        axPV.plot(self['Cavity']['V'][:, 'cRv']*1e6, self['Cavity']['p'][:, 'cRv']/133)
217        axPV.plot(self['Cavity']['V'][:, 'La']*1e6, self['Cavity']['p'][:, 'La']/133)
218        axPV.plot(self['Cavity']['V'][:, 'Ra']*1e6, self['Cavity']['p'][:, 'Ra']/133)
219        axPV.spines[['top', 'right']].set_visible(False)
220        axPV.set_title('Pressure-Volume loop', weight='bold')
221        axPV.set_xlabel('Volume [mL]')
222        axPV.set_ylabel('Pressure [mmHg]')
223        axPV.spines[['bottom', 'left']].set_position(('outward', 5))
224
225        ylabel_x_left = -0.25
226        ylabel_x_right = 1.25
227
228        # Volumes
229        t = self['Solver']['t']*1e3
230
231        axVRv = plt.subplot2grid(grid_size, (0, 0), rowspan=8, colspan=6, fig=fig)
232        axVRv.plot(t, self['Cavity']['V'][:, 'cRv']*1e6)
233        axVRv.plot(t, self['Cavity']['V'][:, 'Ra']*1e6)
234        axVRv.set_ylim(lim_V)
235        axVRv.set_ylabel('Volume\n[mL]')
236        axVRv.spines[['top', 'right']].set_visible(False)
237        axVRv.set_title('Right Heart')
238        # axVRv.set_xticks([])
239        axVRv.tick_params(axis='both', direction='in')
240        axVRv.yaxis.set_label_coords(ylabel_x_left, 0.5)
241
242
243        axVLv = plt.subplot2grid(grid_size, (0, 6), rowspan=8, colspan=6, fig=fig)
244        axVLv.plot(t, self['Cavity']['V'][:, 'cLv']*1e6)
245        axVLv.plot(t, self['Cavity']['V'][:, 'La']*1e6)
246        axVLv.set_ylabel('Volume\n[mL]')
247        axVLv.set_ylim(lim_V)
248        axVLv.yaxis.set_ticks_position('right')
249        axVLv.yaxis.set_label_position('right')
250        axVLv.spines['right'].set_position(('outward', 0))
251        axVLv.spines[['top', 'left']].set_visible(False)
252        axVLv.set_title('Left Heart')
253        # axVLv.set_xticks([])
254        axVLv.tick_params(axis='both', direction='in')
255        axVLv.yaxis.set_label_coords(ylabel_x_right, 0.5)
256
257        # Pressures
258        axpRv = plt.subplot2grid(grid_size, (8, 0), rowspan=8, colspan=6, fig=fig)
259        axpRv.plot(t, self['Cavity']['p'][:, 'cRv']/133)
260        axpRv.plot(t, self['Cavity']['p'][:, 'Ra']/133)
261        axpRv.plot(t, self['Cavity']['p'][:, 'PuArt']/133)
262        axpRv.spines[['top', 'right']].set_visible(False)
263        # axpRv.set_xticks([])
264        axpRv.tick_params(axis='both', direction='in')
265        axpRv.set_ylim(lim_p)
266        axpRv.set_ylabel('Pressure\n[mmHg]')
267        axpRv.yaxis.set_label_coords(ylabel_x_left, 0.5)
268
269        axpLv = plt.subplot2grid(grid_size, (8, 6), rowspan=8, colspan=6, fig=fig)
270        axpLv.plot(t, self['Cavity']['p'][:, 'cLv']/133)
271        axpLv.plot(t, self['Cavity']['p'][:, 'La']/133)
272        axpLv.plot(t, self['Cavity']['p'][:, 'SyArt']/133)
273        axpLv.yaxis.set_ticks_position('right')
274        axpLv.yaxis.set_label_position('right')
275        axpLv.spines['right'].set_position(('outward', 0))
276        axpLv.spines[['top', 'left']].set_visible(False)
277        # axpLv.set_xticks([])
278        axpLv.tick_params(axis='both', direction='in')
279        axpLv.set_ylim(lim_p)
280        axpLv.set_ylabel('Pressure\n[mmHg]')
281        axpLv.yaxis.set_label_coords(ylabel_x_right, 0.5)
282
283        # Valves
284        ax = plt.subplot2grid(grid_size, (16, 0), rowspan=6, colspan=6, fig=fig)
285        ax.plot(t, self['Valve']['q'][:, 'RaRv']*1e6)
286        ax.plot(t, self['Valve']['q'][:, 'RvPuArt']*1e6)
287        ax.spines[['top', 'right']].set_visible(False)
288        ax.set_ylim(lim_q)
289        # ax.set_xticks([])
290        ax.set_ylabel('Flow\n[mL/s]')
291        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
292
293        ax = plt.subplot2grid(grid_size, (16, 6), rowspan=6, colspan=6, fig=fig)
294        ax.plot(t, self['Valve']['q'][:, 'LaLv']*1e6)
295        ax.plot(t, self['Valve']['q'][:, 'LvSyArt']*1e6)
296        ax.spines[['top', 'left']].set_visible(False)
297        ax.set_ylim(lim_q)
298        ax.yaxis.set_ticks_position('right')
299        ax.yaxis.set_label_position('right')
300        # ax.set_xticks([])
301        ax.set_ylabel('Flow\n[mL/s]')
302        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
303
304        # Stress
305        ax = plt.subplot2grid(grid_size, (22, 0), rowspan=4, colspan=6, fig=fig)
306        ax.plot(t, self['Patch']['Sf'][:, 'pRv0']*1e-3)
307        ax.plot(t, self['Patch']['Sf'][:, 'pRa0']*1e-3)
308        ax.spines[['top', 'right']].set_visible(False)
309        # ax.set_xticks([])
310        ax.set_ylim(lim_Sf)
311        ax.set_ylabel('Total\nstress [kPa]')
312        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
313
314        ax = plt.subplot2grid(grid_size, (22, 6), rowspan=4, colspan=6, fig=fig)
315        ax.plot(t, self['Patch']['Sf'][:, 'pLv0']*1e-3)
316        ax.plot(t, self['Patch']['Sf'][:, 'pSv0']*1e-3)
317        ax.plot(t, self['Patch']['Sf'][:, 'pLa0']*1e-3)
318        ax.spines[['top', 'left']].set_visible(False)
319        ax.yaxis.set_ticks_position('right')
320        ax.yaxis.set_label_position('right')
321        # ax.set_xticks([])
322        ax.set_ylim(lim_Sf)
323        ax.set_ylabel('Total\nstress [kPa]')
324        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
325
326        # Sarcomere Length
327        ax = plt.subplot2grid(grid_size, (26, 0), rowspan=6, colspan=6, fig=fig)
328        ax.plot(t, self['Patch']['l_s'][:, 'pRv0'])
329        ax.plot(t, self['Patch']['l_s'][:, 'pRa0'])
330        ax.spines[['top', 'right']].set_visible(False)
331        ax.spines['bottom'].set_position(('outward', 5))
332        ax.set_ylim(lim_Ls)
333        ax.set_ylabel('Sarcomere\nlength [$\mu$m]')
334        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
335
336        ax = plt.subplot2grid(grid_size, (26, 6), rowspan=6, colspan=6, fig=fig)
337        ax.plot(t, self['Patch']['l_s'][:, 'pLv0'])
338        ax.plot(t, self['Patch']['l_s'][:, 'pSv0'])
339        ax.plot(t, self['Patch']['l_s'][:, 'pLa0'])
340        ax.spines[['top', 'left']].set_visible(False)
341        ax.spines['bottom'].set_position(('outward', 5))
342        ax.yaxis.set_ticks_position('right')
343        ax.yaxis.set_label_position('right')
344        ax.set_ylim(lim_Ls)
345        ax.set_ylabel('Sarcomere\nlength [$\mu$m]')
346        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
347
348        # ax.set_xlabel('Time [ms]')
349        # ax.xaxis.set_label_coords(0, -0.3)
350
351
352        # Plot TriSeg
353        titles = ['Pre-A', 'Onset QRS', 'Peak LV \n pressure', 'AV close']
354        idx = [0,
355               np.argmax(np.diff(self['Patch']['C'][:, 'pLv0'])>0),
356               np.argmax(self['Cavity']['p'][:, 'cLv']),
357               len(t) - 1 - np.argmax(
358                   np.diff(self['Valve']['q'][:, 'LvSyArt'][::-1])>0)
359               ]
360
361        for i in range(4):
362            ax = plt.subplot2grid(grid_size, (26, 16+4*i),
363                                  rowspan=5, colspan=4, fig=fig)
364            triseg2022(self, ax, idx[i])
365            ax.spines[['top', 'right', 'bottom', 'left']].set_visible(False)
366            ax.set_xticks([])
367            ax.set_yticks([])
368            plt.xlabel(titles[i], fontsize=12)
369
370        # Plot settings
371        plt.subplots_adjust(
372            top=0.96,
373            bottom=0.05,
374            left=0.075,
375            right=0.98,
376            hspace=5,
377            wspace=0.5)
378        plt.draw()
379
380        # Stress strain
381        ax_right_stress_strain = plt.subplot2grid(grid_size, (18, 17),
382                              rowspan=7, colspan=7, fig=fig)
383        ax_left_stress_strain = plt.subplot2grid(grid_size, (18, 24),
384                              rowspan=7, colspan=7, fig=fig)
385
386        ax_right_stress_strain.plot(
387            self['Patch']['Ef'][:, [1, 4]],
388            self['Patch']['Sf'][:, [1, 4]]*1e-3,
389            )
390        ax_left_stress_strain.plot(
391            self['Patch']['Ef'][:, [0, 2, 3]],
392            self['Patch']['Sf'][:, [0, 2, 3]]*1e-3,
393            )
394
395        ylim = [np.min(self['Patch']['Sf'])*1e-3,
396                np.max(self['Patch']['Sf'])*1e-3]
397        ylim_ptp = np.ptp(ylim)*0.025
398        ylim[0] -= ylim_ptp
399        ylim[1] += ylim_ptp
400        ax_right_stress_strain.set_ylim(ylim)
401        ax_left_stress_strain.set_ylim(ylim)
402        xlim = [np.min(self['Patch']['Ef']),
403                np.max(self['Patch']['Ef'])]
404        xlim_ptp = np.ptp(xlim)*0.025
405        xlim[0] -= xlim_ptp
406        xlim[1] += xlim_ptp
407        ax_right_stress_strain.set_xlim(xlim)
408        ax_left_stress_strain.set_xlim(xlim)
409
410        ax_right_stress_strain.spines[['left', 'bottom']].set_position(('outward', 5))
411        ax_left_stress_strain.spines[['right', 'bottom']].set_position(('outward', 5))
412
413        ax_right_stress_strain.spines[['top', 'right']].set_visible(False)
414        ax_right_stress_strain.set_ylabel('[kPa]', rotation=0, va='bottom', ha='right')
415        ax_right_stress_strain.yaxis.set_label_coords(0., 1.05)
416        ax_right_stress_strain.set_xlabel('Strain [-]')
417        ax_left_stress_strain.spines[['top', 'left']].set_visible(False)
418        ax_left_stress_strain.yaxis.set_ticks_position('right')
419        ax_left_stress_strain.yaxis.set_label_position('right')
420        ax_left_stress_strain.set_ylabel('Stress\n[kPa]', rotation=0, va='bottom', ha='left')
421        ax_left_stress_strain.yaxis.set_label_coords(1., 1.05)
422        ax_left_stress_strain.set_xlabel('Strain [-]')
423
424
425class VanOsta2023_2310(VanOsta2023):
426    def set_reference(self):
427        self['Chamber']['buckling'] = False
428        self['Valve']['soft_closure'] = True
429        self['Valve']['papillary_muscles'] = True
430
431        self['Solver']['dt'] = 0.001
432        self['Solver']['dt_export'] = 0.002
433        self.set('Solver.order',  2)
434
435        self.set('Model.t_cycle', 0.85)
436
437        self.set('Model.PFC.fac',  0.5)
438        self.set('Model.PFC.epsilon', 0.1)
439        self.set('Model.PFC.fac_pfc', 1)
440        self.set('Model.PFC.stable_threshold', 0.001)
441
442        # Set parameters Circulation
443        self['Tube0D']['A0'] = [0.00049828, 0.00049899, 0.00046929, 0.00050574]
444        self['Tube0D']['A_wall'] = [1.13187792e-04, 3.72289316e-05, 7.27489286e-05, 4.31617949e-05]
445        self['Tube0D']['k'] = [1.66666667, 2.33333333, 1.66666667, 2.33333333]
446        self['Tube0D']['l'] = [0.4, 0.4, 0.19, 0.19]
447        self['Tube0D']['p0'] = [12155.25602777,   297.44066943,  1912.6684696 ,   580.51425824]
448
449        self['ArtVen']['p0'] = [6299.86481149, 1000.        ]
450        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
451        self['ArtVen']['k'] = [1.  , 1.88]
452
453        # Set volume state variables
454        # not needed to reset, only if initial model crashes
455        if True:
456            self.set('Model.SyArt.V', 2e-4)
457            self.set('Model.PuArt.V', 1e-4)
458            self.set('Model.SyVen.V', 4e-4)
459            self.set('Model.PuVen.V', 2e-4)
460
461            self.set('Model.Peri.TriSeg.V', 56e-6)
462            self.set('Model.Peri.TriSeg.Y', 36e-3)
463
464            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
465            self.set('Model.Peri.TriSeg.cRv.V',  110e-6)
466            self.set('Model.Peri.La.V',  25e-6)
467            self.set('Model.Peri.Ra.V',  150e-6)
468
469        # all wall parameters
470        self['Patch']['l_se0'] = 0.04
471        self['Patch']['l_s_ref'] = 2.0
472        self['Patch']['l_s0'] = [1.79, 1.79, 1.8 , 1.8 , 1.8 ]
473        self['Patch']['dl_s_pas'] = 0.6
474        self['Patch']['k1'] = 10
475        self['Patch']['dt'] = 0
476        self['Patch']['C_rest'] = 0
477        self['Patch']['l_si0'] = 1.51
478        self['Patch']['LDAD'] = [1.057, 1.057, 0.74 , 0.74 , 0.74 ]
479        self['Patch']['ADO'] = [0.081, 0.081, 0.75 , 0.75 , 0.75 ]
480        self['Patch']['LDCC'] = [4. , 4. , 3., 3., 3.]
481        self['Patch']['v_max'] = 7.
482        self['Patch']['v_max'][:2] = 14.
483
484        self['Patch']['tr'] = [0.4 , 0.4 , 0.21, 0.21, 0.21]
485        self['Patch']['td'] = [0.4 , 0.4 , 0.24, 0.24, 0.24]
486
487        # wall specific
488        self['Patch']['Sf_act'] = [100000., 100000., 120000., 120000., 120000.]
489        self['Patch']['Am_ref'] = [0.00400721, 0.00311263, 0.00918281, 0.00518683, 0.01321999]
490        self['Patch']['V_wall'] = [1.59935455e-05, 7.72140495e-06, 1.16441768e-04, 4.51241560e-05,
491               6.88196851e-05]
492        self['Patch']['Sf_pas'] = [  9.51064014,  10.3647024 , 424.3537027 , 417.74300178,
493               439.21574206]
494
495        # Valves
496        self['Valve']['adaptation_A_open_fac'] = [1.  , 1.12, 1.  , 1.  , 1.12, 1.  ]
497        self['Valve']['A_open'] = [0.00050044, 0.00052598, 0.00046962, 0.00050745, 0.0005583 ,
498               0.00049849]
499        self['Valve']['A_leak'] = [2.64705882e-04, 2.64705882e-10, 2.64705882e-10, 2.64705882e-04,
500               2.64705882e-10, 2.64705882e-10]
501        self['Valve']['l'] = 0.01626978
502        self['Valve']['rho_b'] = 1050
503        self['Valve']['papillary_muscles'] = True
504        self['Valve']['papillary_muscles_slope'] = 100
505        self['Valve']['papillary_muscles_min'] = 0.1
506        self['Valve']['papillary_muscles_A_open_fac'] = 0.1
507        self['Valve']['soft_closure'] = True
508
509        self['Bag']['k'] = 10
510        self['Bag']['V_ref'] = [0.00057937]
511        self['Bag']['p_ref'] = 100
512
513        self.set('Model.Peri.TriSeg.Y', 0.035)
514
515        self['Timings']['law_tau_av'] = 2
516        self['Timings']['c_tau_av0'] = 0.22
517        self['Timings']['c_tau_av1'] = -0.03
518
519        # Adaptation
520        self['Patch']['adapt_gamma'] = 0.5
521        self.set('Solver.store_beats', 1)
522
523        self['PFC']['stable_threshold'] = 1e-3
524
525        self['Patch']['SfPasMaxT'] = [3600., 3600.,  4200.,  4200.,  4200.]
526        self['Patch']['FacSfActT'] = [0.28, 0.28, 0.69, 0.69, 0.69]
527        self['Patch']['SfPasActT'] = [2800., 2800., 6600., 6600., 6600.]
528        self['Patch']['LsPasActT'] = [3.  , 3.  , 2.31, 2.31, 2.31]
529
530        self.run(stable=True)
531
532        return
533
534    def get_unittest_targets(self):
535        """Hardcoded results after initializing and running 1 beat."""
536        return {
537            'LVEDV': 121.75,
538            'LVESV':  48.0,
539            }
540
541
542class VanOsta2023_2306(VanOsta2023):
543    def set_reference(self):
544        self['Chamber']['buckling'] = False
545        self['Valve']['soft_closure'] = True
546
547        self['Solver']['dt'] = 0.0005
548        self['Solver']['dt_export'] = self['Solver']['dt']
549
550        self.set('Model.t_cycle', 0.85)
551
552        self.set('Model.PFC.fac',  0.5)
553        self.set('Model.PFC.epsilon', 0.1)
554        self.set('Model.PFC.fac_pfc', 1)
555        self.set('Model.PFC.stable_threshold', 0.001)
556
557        # Set parameters Circulation
558        self['Tube0D']['A0'] = [
559            0.00049833,
560            0.00050034,
561            0.0004684 ,
562            0.00051849,
563            ]
564        self['Tube0D']['A_wall'] = [
565            1.14408188e-04,
566            4.81082964e-05,
567            5.19472921e-05,
568            5.79933632e-05,
569            ]
570        self['Tube0D']['k'] = np.array([ 8., 10.,  8., 10.]) / 3 - 1
571        self['Tube0D']['l'] = [0.4, 0.4, 0.2, 0.2]
572        self['Tube0D']['p0'] = [
573            12161.21763497,
574            112.09788908,
575            1899.95009955,
576            572.39790382,
577            ]
578
579        self['ArtVen']['p0'] = [6395.43973044,  500.]
580        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
581        self['ArtVen']['k'] = [1, 2]
582
583        # Set volume state variables
584        # not needed to reset, only if initial model crashes
585        if True:
586            self.set('Model.SyArt.V', 2e-4)
587            self.set('Model.PuArt.V', 1e-4)
588            self.set('Model.SyVen.V', 4e-4)
589            self.set('Model.PuVen.V', 2e-4)
590
591            self.set('Model.Peri.TriSeg.V', 44e-6)
592            self.set('Model.Peri.TriSeg.Y', 34.6e-3)
593
594            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
595            self.set('Model.Peri.TriSeg.cRv.V',  100e-6)
596            self.set('Model.Peri.La.V',  100e-6)
597            self.set('Model.Peri.Ra.V',  50e-6)
598
599        # all wall parameters
600        self['Patch']['l_se'] = 0.04
601        self['Patch']['l_s_ref'] = 2.0
602        self['Patch']['l_s0'] = 1.8
603        self['Patch']['dl_s_pas'] = 0.6
604        self['Patch']['k1'] = 10
605        self['Patch']['dt'] = 0
606        self['Patch']['C_rest'] = 0
607        self['Patch']['l_si0'] = 1.51
608        self['Patch']['LDAD'] = 1.057
609        self['Patch']['ADO'] = 0.65
610        self['Patch']['LDCC'] = 4.
611        self['Patch']['v_max'] = 7.
612        self['Patch']['v_max'][:2] = 14.
613
614        self['Patch']['tr'] = 0.25
615        self['Patch']['td'] = 0.25
616        self['Patch']['tr'][:2] = 0.4
617        self['Patch']['td'][:2] = 0.4
618
619        # wall specific
620        self['Patch']['Am_ref'] = [
621            0.00706159,
622            0.00603949,
623            0.0093191 ,
624            0.00507053,
625            0.01313792,
626            ]
627        self['Patch']['V_wall'] = [
628            1.86518494e-05,
629            6.12511827e-06,
630            9.51121820e-05,
631            3.74835076e-05,
632            4.97478018e-05,
633            ]
634        self['Patch']['Sf_pas'] = np.array([
635            51591.64007141,
636            52277.48834349,
637            22370.46712926,
638            22051.80366169,
639            23278.22158555,
640            ]) * 0.0349
641        self['Patch']['Sf_act'] = [84e3, 84e3, 120e3, 120e3, 120e3]
642
643        # Valves
644        self['Valve']['adaptation_A_open_fac'] = [
645            1.,
646            1.5,
647            1.,
648            1.,
649            1.5,
650            1.]
651        self['Valve']['A_open'] = [
652            0.00050078,
653            0.0007021,
654            0.00046806,
655            0.00051776,
656            0.00074755,
657            0.00049837]
658        self['Valve']['A_leak'] = [
659            2.64705882e-04,
660            2.64705882e-10,
661            2.64705882e-10,
662            2.64705882e-04,
663            2.64705882e-10,
664            2.64705882e-10]
665        self['Valve']['l'] = 0.01626978
666        self['Valve']['rho_b'] = 1050
667        self['Valve']['papillary_muscles'] = False
668        self['Valve']['papillary_muscles_slope'] = 100
669        self['Valve']['papillary_muscles_min'] = 0.1
670        self['Valve']['papillary_muscles_A_open_fac'] = 0.1
671        self['Valve']['soft_closure'] = True
672
673        self['Bag']['k'] = 10
674        self['Bag']['V_ref'] = 0.0005
675        self['Bag']['p_ref'] = 100
676
677        self.set('Model.Peri.TriSeg.Y', 0.035)
678
679        self['Timings']['law_tau_av'] = 2
680        self['Timings']['c_tau_av0'] = 0.2
681        self['Timings']['c_tau_av1'] = -0.485*60e-3
682
683        # Adaptation
684        self['Patch']['adapt_gamma'] = 0.75
685
686        self.run(stable=True)
687
688    def get_unittest_targets(self):
689        """Hardcoded results after initializing and running 1 beat."""
690        return {
691            'LVEDV': 141.6,
692            'LVESV':  69.4,
693            }
694
695
696class VanOsta2023_2308(VanOsta2023):
697    def set_reference(self):
698        self['Chamber']['buckling'] = False
699        self['Valve']['soft_closure'] = True
700        self['Valve']['papillary_muscles'] = True
701
702        self['Solver']['dt'] = 0.001
703        self['Solver']['dt_export'] = 0.002
704        self.set('Solver.order',  2)
705
706        self.set('Model.t_cycle', 0.85)
707
708        self.set('Model.PFC.fac',  0.5)
709        self.set('Model.PFC.epsilon', 0.1)
710        self.set('Model.PFC.fac_pfc', 1)
711        self.set('Model.PFC.stable_threshold', 0.001)
712
713        # Set parameters Circulation
714        self['Tube0D']['A0'] = [0.0004982 , 0.00049959, 0.00045184, 0.00051768]
715        self['Tube0D']['A_wall'] = [1.13597208e-04, 3.79106863e-05, 8.91394588e-05, 4.26698697e-05]
716        self['Tube0D']['k'] = [1.66666667, 2.33333333, 1.        , 2.33333333]
717        self['Tube0D']['l'] = [0.4, 0.4, 0.1, 0.1]
718        self['Tube0D']['p0'] = [12154.79898845,   213.26733196,  1913.05480512,   600.72294959]
719
720        self['ArtVen']['p0'] = [6345.26731406,  950.        ]
721        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
722        self['ArtVen']['k'] = [1.  , 1.72]
723
724        # Set volume state variables
725        # not needed to reset, only if initial model crashes
726        if True:
727            self.set('Model.SyArt.V', 2e-4)
728            self.set('Model.PuArt.V', 1e-4)
729            self.set('Model.SyVen.V', 4e-4)
730            self.set('Model.PuVen.V', 2e-4)
731
732            self.set('Model.Peri.TriSeg.V', 56e-6)
733            self.set('Model.Peri.TriSeg.Y', 36e-3)
734
735            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
736            self.set('Model.Peri.TriSeg.cRv.V',  110e-6)
737            self.set('Model.Peri.La.V',  25e-6)
738            self.set('Model.Peri.Ra.V',  150e-6)
739
740        # all wall parameters
741        self['Patch']['l_se'] = 0.04
742        self['Patch']['l_s_ref'] = 2.0
743        self['Patch']['l_s0'] = 1.8
744        self['Patch']['dl_s_pas'] = 0.6
745        self['Patch']['k1'] = 10
746        self['Patch']['dt'] = 0
747        self['Patch']['C_rest'] = 0
748        self['Patch']['l_si0'] = 1.51
749        self['Patch']['LDAD'] = [1.057, 1.057, 0.64 , 0.64 , 0.64 ]
750        self['Patch']['ADO'] = [0.65, 0.65, 0.75, 0.75, 0.75]
751        self['Patch']['LDCC'] = [4. , 4. , 3.2, 3.2, 3.2]
752        self['Patch']['v_max'] = 7.
753        self['Patch']['v_max'][:2] = 14.
754
755        self['Patch']['tr'] = [0.4 , 0.4 , 0.24, 0.24, 0.24]
756        self['Patch']['td'] = [0.4 , 0.4 , 0.23, 0.23, 0.23]
757
758        # wall specific
759        self['Patch']['Sf_act'] = [80e3, 80e3, 120e3, 120e3, 120e3]
760        self['Patch']['Am_ref'] = [0.00458522, 0.00381765, 0.00806997, 0.00454307, 0.01176514]
761        self['Patch']['V_wall'] = [2.75050615e-05, 1.25959326e-05, 8.80910385e-05, 3.46281364e-05,
762               5.16086110e-05]
763        self['Patch']['Sf_pas'] = [ 16.51117993,  17.0639141 , 580.39235927, 569.67392472,
764               606.50514861]
765
766        # Valves
767        self['Valve']['adaptation_A_open_fac'] = [1.  , 1.11, 1.  , 1.  , 1.11, 1.  ]
768        self['Valve']['A_open'] = [0.00050042, 0.00050203, 0.00045228, 0.00051764, 0.00055321,
769               0.00049838]
770        self['Valve']['A_leak'] = [2.64705882e-04, 2.64705882e-10, 2.64705882e-10, 2.64705882e-04,
771               2.64705882e-10, 2.64705882e-10]
772        self['Valve']['l'] = 0.01626978
773        self['Valve']['rho_b'] = 1050
774        self['Valve']['papillary_muscles'] = True
775        self['Valve']['papillary_muscles_slope'] = 100
776        self['Valve']['papillary_muscles_min'] = 0.1
777        self['Valve']['papillary_muscles_A_open_fac'] = 0.1
778        self['Valve']['soft_closure'] = True
779
780        self['Bag']['k'] = 10
781        self['Bag']['V_ref'] = [0.00051968]
782        self['Bag']['p_ref'] = 100
783
784        self.set('Model.Peri.TriSeg.Y', 0.035)
785
786        self['Timings']['law_tau_av'] = 2
787        self['Timings']['c_tau_av0'] = 0.172
788        self['Timings']['c_tau_av1'] = -0.485*60e-3
789
790        # Adaptation
791        self['Patch']['adapt_gamma'] = 0.5
792        self.set('Solver.store_beats', 1)
793
794        self['PressureFlowControl']['stable_threshold'] = 1e-3
795
796        self['Patch']['SfPasMaxT'] = [6400., 6400.,  6600.,  6600.,  6600.]
797        self['Patch']['FacSfActT'] = [0.44, 0.44, 0.61, 0.61, 0.61]
798        self['Patch']['SfPasActT'] = [4800., 4800., 6600., 6600., 6600.]
799        self['Patch']['LsPasActT'] = [3.  , 3.  , 2.23, 2.23, 2.23]
800
801        self.run(stable=True)
802
803        return
804
805    def get_unittest_targets(self):
806        """Hardcoded results after initializing and running 1 beat."""
807        return {
808            'LVEDV': 119.7,
809            'LVESV':  47.2,
810            }
811
812
813class Walmsley2015(VanOsta2023):
814    def set_reference(self):
815
816        # change model parameters here
817        self['Chamber']['buckling'] = False
818        self['Valve']['soft_closure'] = True
819        self['Valve']['papillary_muscles'] = False
820
821        self['Solver']['dt'] = 0.001
822        self['Solver']['dt_export'] = 0.002
823        self.set('Solver.order',  2)
824
825        self.set('Model.t_cycle', 0.85)
826
827        self.set('Model.PFC.fac',  0.5)
828        self.set('Model.PFC.epsilon', 0.1)
829        self.set('Model.PFC.fac_pfc', 1)
830        self.set('Model.PFC.stable_threshold', 0.001)
831
832        # Set parameters Circulation
833        self['Tube0D']['A0'] = [4.9720e-04, 4.9830e-04, 4.6850e-04, 5.0910e-04]
834        self['Tube0D']['A_wall'] = [1.1370e-04, 6.410e-05, 9.670e-05, 5.950e-05]
835        self['Tube0D']['k'] = [1.66666667, 2.33333333, 1.66666667, 2.33333333]
836        # self['Tube0D']['k'] = [8, 10, 8, 10]
837        self['Tube0D']['l'] = [0.4, 0.4, 0.2, 0.2]
838        self['Tube0D']['p0'] = [12164,   228,  2630 , 678]
839
840        self['ArtVen']['p0'] = [6396, 1500.        ]
841        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
842        self['ArtVen']['k'] = [1.  , 2]
843
844        # Set volume state variables
845        # not needed to reset, only if initial model crashes
846        if True:
847            self.set('Model.SyArt.V', 2e-4)
848            self.set('Model.PuArt.V', 1e-4)
849            self.set('Model.SyVen.V', 4e-4)
850            self.set('Model.PuVen.V', 2e-4)
851
852            self.set('Model.Peri.TriSeg.V', 56e-6)
853            self.set('Model.Peri.TriSeg.Y', 36e-3)
854
855            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
856            self.set('Model.Peri.TriSeg.cRv.V',  110e-6)
857            self.set('Model.Peri.La.V',  25e-6)
858            self.set('Model.Peri.Ra.V',  150e-6)
859
860        # all wall parameters
861        self['Patch']['l_se0'] = 0.04
862        self['Patch']['l_s_ref'] = 2.0
863        self['Patch']['l_s0'] = [1.8, 1.8, 1.8 , 1.8 , 1.8 ]
864        self['Patch']['dl_s_pas'] = 0.6
865        self['Patch']['k1'] = 10
866        self['Patch']['dt'] = 0
867        self['Patch']['C_rest'] = 0
868        self['Patch']['l_si0'] = 1.51
869        self['Patch']['LDAD'] = [1.057, 1.057, 0.74 , 0.74 , 0.74 ]
870        self['Patch']['ADO'] = [0.081, 0.081, 0.75 , 0.75 , 0.75 ]
871        self['Patch']['LDCC'] = [4. , 4. , 3., 3., 3.]
872        self['Patch']['v_max'] = 7.
873        self['Patch']['v_max'][:2] = 14.
874
875        self['Patch']['tr'] = [0.4 , 0.4 , 0.25, 0.25, 0.25]
876        self['Patch']['td'] = [0.4 , 0.4 , 0.25, 0.25, 0.25]
877
878        # wall specific
879        self['Patch']['Sf_act'] = [60000., 60000., 100000., 100000., 100000.]
880        self['Patch']['Am_ref'] = [0.0071, 0.0058, 0.0099, 0.0052, 0.0135]
881        self['Patch']['V_wall'] = [1.3377e-05, 5.4588e-06, 1.1382e-04, 3.8320e-05,
882               8.0093e-05]
883        self['Patch']['Sf_pas'] = [  2.8942e3,  3.6579e3 , 794.1495 , 790.6246,
884               801.1993]
885
886        # Valves
887        self['Valve']['adaptation_A_open_fac'] = [1.  , 1.12, 1.  , 1.  , 1.12, 1.  ]
888        self['Valve']['A_open'] = [0.00049811, 0.00070281, 0.00046854, 0.00050903, 0.00074578 ,
889               0.00049719]
890        self['Valve']['A_leak'] = [2.6471e-04, 2.6471e-10, 2.6471e-10, 2.6471e-04,
891               2.6471e-10, 2.6471e-10]
892        self['Valve']['l'] = 0.0163
893        self['Valve']['rho_b'] = 1050
894        self['Valve']['papillary_muscles'] = False
895        self['Valve']['papillary_muscles_slope'] = 100
896        self['Valve']['papillary_muscles_min'] = 0.1
897        self['Valve']['papillary_muscles_A_open_fac'] = 0.1
898        self['Valve']['soft_closure'] = True
899
900        self['Bag']['k'] = 10
901        self['Bag']['V_ref'] = [0.00072379]
902        self['Bag']['p_ref'] = 100
903
904        self.set('Model.Peri.TriSeg.Y', 0.035)
905
906        # self['Timings']['law_tauAv'] = 1
907        # self['Timings']['c_tauAv0'] = 0.22
908        # self['Timings']['c_tauAv1'] = -0.03
909
910        # Adaptation
911        # self['Patch2022']['adapt_gamma'] = 0.5
912        self.set('Solver.store_beats', 1)
913
914        self['PFC']['stable_threshold'] = 1e-3
915
916        # self['Patch']['SfPasMaxT'] = [3600., 3600.,  4200.,  4200.,  4200.]
917        # self['Patch']['FacSfActT'] = [0.28, 0.28, 0.69, 0.69, 0.69]
918        # self['Patch']['SfPasActT'] = [2800., 2800., 6600., 6600., 6600.]
919        # self['Patch']['LsPasActT'] = [3.  , 3.  , 2.31, 2.31, 2.31]
920
921        self.run(stable=True)