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
  6import circadapt
  7import circadapt.plot
  8from circadapt.model import Model
  9from circadapt.adapt import ModelAdapt
 10from circadapt.plot import triseg2022, mmode
 11import matplotlib.pyplot as plt
 12import numpy as np
 13
 14class VanOsta2023(Model, ModelAdapt):
 15    def __init__(self,
 16                 solver: str = None,
 17                 path_to_circadapt: str = None,
 18                 model_state: dict = None,
 19                 ):
 20        if solver is None:
 21            solver = 'backward_differential'
 22
 23        self._local_save_reference = True
 24
 25        ModelAdapt.__init__(self)
 26        Model.__init__(self,
 27                       solver,
 28                       path_to_circadapt=path_to_circadapt,
 29                       model_state=model_state,
 30                       )
 31
 32    def build(self):
 33        # Circulation
 34        self.add_smart_component('ArtVen', build='SystemicCirculation')
 35        self.add_smart_component('ArtVen', build='PulmonaryCirculation')
 36        self.add_smart_component('Heart', patch_type='Patch2022')
 37        self.add_smart_component('Timings')
 38        self.add_smart_component('PressureFlowControl')
 39
 40        # manually set papillary muscles
 41        self.set_component("Peri.RaRv.wPapMus", "Peri.TriSeg.wRv")
 42        self.set_component("Peri.LaLv.wPapMus", "Peri.TriSeg.wLv")
 43
 44    def set_reference(self):
 45        self['Chamber2022']['buckling'] = False
 46        self['Valve2022']['soft_closure'] = True
 47        self['Valve2022']['papillary_muscles'] = True
 48
 49        self['Solver']['dt'] = 0.001
 50        self['Solver']['dt_export'] = 0.002
 51        self.set('Solver.order',  2)
 52
 53        self.set('Model.t_cycle', 0.85)
 54
 55        self.set('Model.PFC.fac',  0.5)
 56        self.set('Model.PFC.epsilon', 0.1)
 57        self.set('Model.PFC.fac_pfc', 1)
 58        self.set('Model.PFC.stable_threshold', 0.001)
 59
 60        # Set parameters Circulation
 61        self['Tube0D']['A0'] = [0.00049828, 0.00049899, 0.00046929, 0.00050574]
 62        self['Tube0D']['A_wall'] = [1.13187792e-04, 3.72289316e-05, 7.27489286e-05, 4.31617949e-05]
 63        self['Tube0D']['k'] = [1.66666667, 2.33333333, 1.66666667, 2.33333333]
 64        self['Tube0D']['l'] = [0.4, 0.4, 0.19, 0.19]
 65        self['Tube0D']['p0'] = [12155.25602777,   297.44066943,  1912.6684696 ,   580.51425824]
 66
 67        self['ArtVen']['p0'] = [6299.86481149, 1000.        ]
 68        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
 69        self['ArtVen']['k'] = [1.  , 1.88]
 70
 71        # Set volume state variables
 72        # not needed to reset, only if initial model crashes
 73        if True:
 74            self.set('Model.SyArt.V', 2e-4)
 75            self.set('Model.PuArt.V', 1e-4)
 76            self.set('Model.SyVen.V', 4e-4)
 77            self.set('Model.PuVen.V', 2e-4)
 78
 79            self.set('Model.Peri.TriSeg.V', 56e-6)
 80            self.set('Model.Peri.TriSeg.Y', 36e-3)
 81
 82            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
 83            self.set('Model.Peri.TriSeg.cRv.V',  110e-6)
 84            self.set('Model.Peri.La.V',  25e-6)
 85            self.set('Model.Peri.Ra.V',  150e-6)
 86
 87        # all wall parameters
 88        self['Patch2022']['l_se0'] = 0.04
 89        self['Patch2022']['l_s_ref'] = 2.0
 90        self['Patch2022']['l_s0'] = [1.79, 1.79, 1.8 , 1.8 , 1.8 ]
 91        self['Patch2022']['dl_s_pas'] = 0.6
 92        self['Patch2022']['k1'] = 10
 93        self['Patch2022']['dt'] = 0
 94        self['Patch2022']['C_rest'] = 0
 95        self['Patch2022']['l_si0'] = 1.51
 96        self['Patch2022']['LDAD'] = [1.057, 1.057, 0.74 , 0.74 , 0.74 ]
 97        self['Patch2022']['ADO'] = [0.081, 0.081, 0.75 , 0.75 , 0.75 ]
 98        self['Patch2022']['LDCC'] = [4. , 4. , 3., 3., 3.]
 99        self['Patch2022']['v_max'] = 7.
100        self['Patch2022']['v_max'][:2] = 14.
101
102        self['Patch2022']['tr'] = [0.4 , 0.4 , 0.21, 0.21, 0.21]
103        self['Patch2022']['td'] = [0.4 , 0.4 , 0.24, 0.24, 0.24]
104
105        # wall specific
106        self['Patch2022']['Sf_act'] = [100000., 100000., 120000., 120000., 120000.]
107        self['Patch2022']['Am_ref'] = [0.00400721, 0.00311263, 0.00918281, 0.00518683, 0.01321999]
108        self['Patch2022']['V_wall'] = [1.59935455e-05, 7.72140495e-06, 1.16441768e-04, 4.51241560e-05,
109               6.88196851e-05]
110        self['Patch2022']['Sf_pas'] = [  9.51064014,  10.3647024 , 424.3537027 , 417.74300178,
111               439.21574206]
112
113        # Valves
114        self['Valve2022']['adaptation_A_open_fac'] = [1.  , 1.12, 1.  , 1.  , 1.12, 1.  ]
115        self['Valve2022']['A_open'] = [0.00050044, 0.00052598, 0.00046962, 0.00050745, 0.0005583 ,
116               0.00049849]
117        self['Valve2022']['A_leak'] = [2.64705882e-04, 2.64705882e-10, 2.64705882e-10, 2.64705882e-04,
118               2.64705882e-10, 2.64705882e-10]
119        self['Valve2022']['l'] = 0.01626978
120        self['Valve2022']['rho_b'] = 1050
121        self['Valve2022']['papillary_muscles'] = True
122        self['Valve2022']['papillary_muscles_slope'] = 100
123        self['Valve2022']['papillary_muscles_min'] = 0.1
124        self['Valve2022']['papillary_muscles_A_open_fac'] = 0.1
125        self['Valve2022']['soft_closure'] = True
126
127        self['Bag']['k'] = 10
128        self['Bag']['V_ref'] = [0.00057937]
129        self['Bag']['p_ref'] = 100
130
131        self.set('Model.Peri.TriSeg.Y', 0.035)
132
133        self['Timings']['law_tauAv'] = 2
134        self['Timings']['c_tauAv0'] = 0.22
135        self['Timings']['c_tauAv1'] = -0.03
136
137        # Adaptation
138        self['Patch2022']['adapt_gamma'] = 0.5
139        self.set('Solver.store_beats', 1)
140
141        self['PressureFlowControl']['stable_threshold'] = 1e-3
142
143        self['Patch2022']['SfPasMaxT'] = [3600., 3600.,  4200.,  4200.,  4200.]
144        self['Patch2022']['FacSfActT'] = [0.28, 0.28, 0.69, 0.69, 0.69]
145        self['Patch2022']['SfPasActT'] = [2800., 2800., 6600., 6600., 6600.]
146        self['Patch2022']['LsPasActT'] = [3.  , 3.  , 2.31, 2.31, 2.31]
147
148        self.run(stable=True)
149
150        return
151
152
153        options = self.get_adapt_options()
154        self['General']['q0'] = options['exercise']['q0']
155        self['General']['t_cycle'] = options['exercise']['t_cycle']
156        self.run(stable=True)
157
158        reference_y = np.array([
159            self['Patch2022']['SfEcmMax'],
160            self['Patch2022']['SfActMax'],
161            self['Patch2022']['SfPasAct'],
162            self['Patch2022']['LsPasAct'],
163            ]).T
164
165        targets = np.array([
166            np.mean(reference_y[:2, :], axis=0),
167            np.mean(reference_y[2:, :], axis=0),
168            ])
169        targets[:, 1] /= self['Patch2022']['Sf_act'][:][[0, 2]]
170
171        # targets = np.array(
172        #     [[5.5e+05, 1.9e-01, 3.5e+03, 2.2e+00],
173        #      [6.0e+03, 4.9e-01, 4.0e+03, 2.3e+00],
174        #      ])
175
176        self['Patch2022']['SfPasMaxT'] = targets[[0, 0, 1, 1, 1], 0]
177        self['Patch2022']['FacSfActT'] = targets[[0, 0, 1, 1, 1], 1]
178        self['Patch2022']['SfPasActT'] = targets[[0, 0, 1, 1, 1], 2]
179        self['Patch2022']['LsPasActT'] = targets[[0, 0, 1, 1, 1], 3]
180
181        self['Patch2022']['adapt_gamma'] = 0.1
182        self.adapt(verbose=True)
183        self['Patch2022']['adapt_gamma'] = 0.5
184        self.adapt(verbose=True)
185
186        # Set adaptation constants
187        self.calculate_and_set_matrix(verbose=True)
188
189        # set resting state
190        self['General']['q0'] = options['rest']['q0']
191        self['General']['t_cycle'] = options['rest']['t_cycle']
192        self.run(stable=True)
193
194    def get_unittest_targets(self):
195        """Hardcoded results after initializing and running 1 beat."""
196        return {
197            'LVEDV': 121.75,
198            'LVESV':  48.0,
199            }
200
201    def get_unittest_results(self, model):
202        """Real-time results after initializing and running 1 beat."""
203        LVEDV = np.max(model['Cavity']['V'][:, 'cLv'])*1e6
204        LVESV = np.min(model['Cavity']['V'][:, 'cLv'])*1e6
205        return {
206            'LVEDV': LVEDV,
207            'LVESV': LVESV,
208            }
209
210    def plot(self, fig=None):
211        # TODO
212        self.plot_extended(fig)
213
214    def plot_extended(self, fig=None):
215        if fig is None:
216            fig = 1
217        if isinstance(fig, int):
218            fig = plt.figure(fig, clear=True, figsize=(12, 8))
219
220        # Settings
221        grid_size = [32, 32]
222
223        def get_lim(module, signal, locs=slice(None, None, None)):
224            signal = self[module][signal][:, locs]
225            lim = np.array([np.min(signal), np.max(signal)])
226            lim += np.array([-1, 1]) * 0.1*np.diff(lim)
227            return lim
228
229        lim_V = get_lim('Cavity', 'V', ['cRv', 'Ra', 'La', 'cLv']) * 1e6
230        lim_V[0] = 0
231        lim_p = get_lim('Cavity', 'p', ['cRv', 'Ra', 'La', 'cLv']) / 133
232        lim_p[0] = np.min([lim_p[0], 0])
233        lim_Ls = get_lim('Patch2022', 'l_s')
234        lim_Sf = get_lim('Patch2022', 'Sf') * 1e-3
235        lim_q = get_lim('Valve2022', 'q',
236                        ['LaLv', 'RaRv', 'LvSyArt', 'RaPuArt']) * 1e6
237
238        all_lim = [lim_V, lim_p, lim_Ls, lim_Sf, lim_q]
239        if (np.any(np.isnan(all_lim)) or np.any(np.isinf(all_lim))):
240            lim_V = [0, 200]
241            lim_p = [0, 150]
242            lim_Ls = [1.5, 2.0]
243            lim_Sf = [0, 100]
244            lim_q = [-1e-3, 1e-3]
245
246        # Pressure Volume plot
247        axPV = plt.subplot2grid(grid_size, (0, 17), rowspan=15, colspan=15, fig=fig)
248        axPV.plot(self['Cavity']['V'][:, 'cLv']*1e6, self['Cavity']['p'][:, 'cLv']/133)
249        axPV.plot(self['Cavity']['V'][:, 'cRv']*1e6, self['Cavity']['p'][:, 'cRv']/133)
250        axPV.plot(self['Cavity']['V'][:, 'La']*1e6, self['Cavity']['p'][:, 'La']/133)
251        axPV.plot(self['Cavity']['V'][:, 'Ra']*1e6, self['Cavity']['p'][:, 'Ra']/133)
252        axPV.spines[['top', 'right']].set_visible(False)
253        axPV.set_title('Pressure-Volume loop', weight='bold')
254        axPV.set_xlabel('Volume [mL]')
255        axPV.set_ylabel('Pressure [mmHg]')
256        axPV.spines[['bottom', 'left']].set_position(('outward', 5))
257
258        ylabel_x_left = -0.25
259        ylabel_x_right = 1.25
260
261        # Volumes
262        t = self['Solver']['t']*1e3
263
264        axVRv = plt.subplot2grid(grid_size, (0, 0), rowspan=8, colspan=6, fig=fig)
265        axVRv.plot(t, self['Cavity']['V'][:, 'cRv']*1e6)
266        axVRv.plot(t, self['Cavity']['V'][:, 'Ra']*1e6)
267        axVRv.set_ylim(lim_V)
268        axVRv.set_ylabel('Volume\n[mL]')
269        axVRv.spines[['top', 'right']].set_visible(False)
270        axVRv.set_title('Right Heart')
271        # axVRv.set_xticks([])
272        axVRv.tick_params(axis='both', direction='in')
273        axVRv.yaxis.set_label_coords(ylabel_x_left, 0.5)
274
275
276        axVLv = plt.subplot2grid(grid_size, (0, 6), rowspan=8, colspan=6, fig=fig)
277        axVLv.plot(t, self['Cavity']['V'][:, 'cLv']*1e6)
278        axVLv.plot(t, self['Cavity']['V'][:, 'La']*1e6)
279        axVLv.set_ylabel('Volume\n[mL]')
280        axVLv.set_ylim(lim_V)
281        axVLv.yaxis.set_ticks_position('right')
282        axVLv.yaxis.set_label_position('right')
283        axVLv.spines['right'].set_position(('outward', 0))
284        axVLv.spines[['top', 'left']].set_visible(False)
285        axVLv.set_title('Left Heart')
286        # axVLv.set_xticks([])
287        axVLv.tick_params(axis='both', direction='in')
288        axVLv.yaxis.set_label_coords(ylabel_x_right, 0.5)
289
290        # Pressures
291        axpRv = plt.subplot2grid(grid_size, (8, 0), rowspan=8, colspan=6, fig=fig)
292        axpRv.plot(t, self['Cavity']['p'][:, 'cRv']/133)
293        axpRv.plot(t, self['Cavity']['p'][:, 'Ra']/133)
294        axpRv.plot(t, self['Cavity']['p'][:, 'PuArt']/133)
295        axpRv.spines[['top', 'right']].set_visible(False)
296        # axpRv.set_xticks([])
297        axpRv.tick_params(axis='both', direction='in')
298        axpRv.set_ylim(lim_p)
299        axpRv.set_ylabel('Pressure\n[mmHg]')
300        axpRv.yaxis.set_label_coords(ylabel_x_left, 0.5)
301
302        axpLv = plt.subplot2grid(grid_size, (8, 6), rowspan=8, colspan=6, fig=fig)
303        axpLv.plot(t, self['Cavity']['p'][:, 'cLv']/133)
304        axpLv.plot(t, self['Cavity']['p'][:, 'La']/133)
305        axpLv.plot(t, self['Cavity']['p'][:, 'SyArt']/133)
306        axpLv.yaxis.set_ticks_position('right')
307        axpLv.yaxis.set_label_position('right')
308        axpLv.spines['right'].set_position(('outward', 0))
309        axpLv.spines[['top', 'left']].set_visible(False)
310        # axpLv.set_xticks([])
311        axpLv.tick_params(axis='both', direction='in')
312        axpLv.set_ylim(lim_p)
313        axpLv.set_ylabel('Pressure\n[mmHg]')
314        axpLv.yaxis.set_label_coords(ylabel_x_right, 0.5)
315
316        # Valves
317        ax = plt.subplot2grid(grid_size, (16, 0), rowspan=6, colspan=6, fig=fig)
318        ax.plot(t, self['Valve2022']['q'][:, 'RaRv']*1e6)
319        ax.plot(t, self['Valve2022']['q'][:, 'RvPuArt']*1e6)
320        ax.spines[['top', 'right']].set_visible(False)
321        ax.set_ylim(lim_q)
322        # ax.set_xticks([])
323        ax.set_ylabel('Flow\n[mL/s]')
324        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
325
326        ax = plt.subplot2grid(grid_size, (16, 6), rowspan=6, colspan=6, fig=fig)
327        ax.plot(t, self['Valve2022']['q'][:, 'LaLv']*1e6)
328        ax.plot(t, self['Valve2022']['q'][:, 'LvSyArt']*1e6)
329        ax.spines[['top', 'left']].set_visible(False)
330        ax.set_ylim(lim_q)
331        ax.yaxis.set_ticks_position('right')
332        ax.yaxis.set_label_position('right')
333        # ax.set_xticks([])
334        ax.set_ylabel('Flow\n[mL/s]')
335        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
336
337        # Stress
338        ax = plt.subplot2grid(grid_size, (22, 0), rowspan=4, colspan=6, fig=fig)
339        ax.plot(t, self['Patch2022']['Sf'][:, 'pRv0']*1e-3)
340        ax.plot(t, self['Patch2022']['Sf'][:, 'pRa0']*1e-3)
341        ax.spines[['top', 'right']].set_visible(False)
342        # ax.set_xticks([])
343        ax.set_ylim(lim_Sf)
344        ax.set_ylabel('Total\nstress [kPa]')
345        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
346
347        ax = plt.subplot2grid(grid_size, (22, 6), rowspan=4, colspan=6, fig=fig)
348        ax.plot(t, self['Patch2022']['Sf'][:, 'pLv0']*1e-3)
349        ax.plot(t, self['Patch2022']['Sf'][:, 'pSv0']*1e-3)
350        ax.plot(t, self['Patch2022']['Sf'][:, 'pLa0']*1e-3)
351        ax.spines[['top', 'left']].set_visible(False)
352        ax.yaxis.set_ticks_position('right')
353        ax.yaxis.set_label_position('right')
354        # ax.set_xticks([])
355        ax.set_ylim(lim_Sf)
356        ax.set_ylabel('Total\nstress [kPa]')
357        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
358
359        # Sarcomere Length
360        ax = plt.subplot2grid(grid_size, (26, 0), rowspan=6, colspan=6, fig=fig)
361        ax.plot(t, self['Patch2022']['l_s'][:, 'pRv0'])
362        ax.plot(t, self['Patch2022']['l_s'][:, 'pRa0'])
363        ax.spines[['top', 'right']].set_visible(False)
364        ax.spines['bottom'].set_position(('outward', 5))
365        ax.set_ylim(lim_Ls)
366        ax.set_ylabel('Sarcomere\nlength [$\mu$m]')
367        ax.yaxis.set_label_coords(ylabel_x_left, 0.5)
368
369        ax = plt.subplot2grid(grid_size, (26, 6), rowspan=6, colspan=6, fig=fig)
370        ax.plot(t, self['Patch2022']['l_s'][:, 'pLv0'])
371        ax.plot(t, self['Patch2022']['l_s'][:, 'pSv0'])
372        ax.plot(t, self['Patch2022']['l_s'][:, 'pLa0'])
373        ax.spines[['top', 'left']].set_visible(False)
374        ax.spines['bottom'].set_position(('outward', 5))
375        ax.yaxis.set_ticks_position('right')
376        ax.yaxis.set_label_position('right')
377        ax.set_ylim(lim_Ls)
378        ax.set_ylabel('Sarcomere\nlength [$\mu$m]')
379        ax.yaxis.set_label_coords(ylabel_x_right, 0.5)
380
381        # ax.set_xlabel('Time [ms]')
382        # ax.xaxis.set_label_coords(0, -0.3)
383
384
385        # Plot TriSeg
386        titles = ['Pre-A', 'Onset QRS', 'Peak LV \n pressure', 'AV close']
387        idx = [0,
388               np.argmax(np.diff(self['Patch2022']['C'][:, 'pLv0'])>0),
389               np.argmax(self['Cavity']['p'][:, 'cLv']),
390               len(t) - 1 - np.argmax(
391                   np.diff(self['Valve2022']['q'][:, 'LvSyArt'][::-1])>0)
392               ]
393
394        for i in range(4):
395            ax = plt.subplot2grid(grid_size, (26, 16+4*i),
396                                  rowspan=5, colspan=4, fig=fig)
397            triseg2022(self, ax, idx[i])
398            ax.spines[['top', 'right', 'bottom', 'left']].set_visible(False)
399            ax.set_xticks([])
400            ax.set_yticks([])
401            plt.xlabel(titles[i], fontsize=12)
402
403        # Plot settings
404        plt.subplots_adjust(
405            top=0.96,
406            bottom=0.05,
407            left=0.075,
408            right=0.98,
409            hspace=5,
410            wspace=0.5)
411        plt.draw()
412
413        # Plot MMode
414        ax_mmode = plt.subplot2grid(grid_size, (17, 17),
415                              rowspan=8, colspan=15, fig=fig)
416        circadapt.plot.triseg.mmode(self, ax_mmode)
417        ax_mmode.axhline(0, c='k', ls='--')
418        ax_mmode.spines[['top', 'right']].set_visible(False)
419
420        # Plot Y
421        # ax_mmode = plt.subplot2grid(grid_size, (17, 25),
422        #                       rowspan=8, colspan=7, fig=fig)
423        ax_mmode.plot(t, self['TriSeg2022']['Y']*1e3 * np.array([[1, -1]]), c='k')
424        ax_mmode.spines[['top', 'right']].set_visible(False)
425        plt.ylabel('MMode and Y [mm]')
426
427
428class VanOsta2023_2306(VanOsta2023):
429    def set_reference(self):
430        self['Chamber2022']['buckling'] = False
431        self['Valve2022']['soft_closure'] = True
432
433        self['Solver']['dt'] = 0.0005
434        self['Solver']['dt_export'] = self['Solver']['dt']
435
436        self.set('Model.t_cycle', 0.85)
437
438        self.set('Model.PFC.fac',  0.5)
439        self.set('Model.PFC.epsilon', 0.1)
440        self.set('Model.PFC.fac_pfc', 1)
441        self.set('Model.PFC.stable_threshold', 0.001)
442
443        # Set parameters Circulation
444        self['Tube0D']['A0'] = [
445            0.00049833,
446            0.00050034,
447            0.0004684 ,
448            0.00051849,
449            ]
450        self['Tube0D']['A_wall'] = [
451            1.14408188e-04,
452            4.81082964e-05,
453            5.19472921e-05,
454            5.79933632e-05,
455            ]
456        self['Tube0D']['k'] = np.array([ 8., 10.,  8., 10.]) / 3 - 1
457        self['Tube0D']['l'] = [0.4, 0.4, 0.2, 0.2]
458        self['Tube0D']['p0'] = [
459            12161.21763497,
460            112.09788908,
461            1899.95009955,
462            572.39790382,
463            ]
464
465        self['ArtVen']['p0'] = [6395.43973044,  500.]
466        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
467        self['ArtVen']['k'] = [1, 2]
468
469        # Set volume state variables
470        # not needed to reset, only if initial model crashes
471        if True:
472            self.set('Model.SyArt.V', 2e-4)
473            self.set('Model.PuArt.V', 1e-4)
474            self.set('Model.SyVen.V', 4e-4)
475            self.set('Model.PuVen.V', 2e-4)
476
477            self.set('Model.Peri.TriSeg.V', 44e-6)
478            self.set('Model.Peri.TriSeg.Y', 34.6e-3)
479
480            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
481            self.set('Model.Peri.TriSeg.cRv.V',  100e-6)
482            self.set('Model.Peri.La.V',  100e-6)
483            self.set('Model.Peri.Ra.V',  50e-6)
484
485        # all wall parameters
486        self['Patch2022']['l_se'] = 0.04
487        self['Patch2022']['l_s_ref'] = 2.0
488        self['Patch2022']['l_s0'] = 1.8
489        self['Patch2022']['dl_s_pas'] = 0.6
490        self['Patch2022']['k1'] = 10
491        self['Patch2022']['dt'] = 0
492        self['Patch2022']['C_rest'] = 0
493        self['Patch2022']['l_si0'] = 1.51
494        self['Patch2022']['LDAD'] = 1.057
495        self['Patch2022']['ADO'] = 0.65
496        self['Patch2022']['LDCC'] = 4.
497        self['Patch2022']['v_max'] = 7.
498        self['Patch2022']['v_max'][:2] = 14.
499
500        self['Patch2022']['tr'] = 0.25
501        self['Patch2022']['td'] = 0.25
502        self['Patch2022']['tr'][:2] = 0.4
503        self['Patch2022']['td'][:2] = 0.4
504
505        # wall specific
506        self['Patch2022']['Am_ref'] = [
507            0.00706159,
508            0.00603949,
509            0.0093191 ,
510            0.00507053,
511            0.01313792,
512            ]
513        self['Patch2022']['V_wall'] = [
514            1.86518494e-05,
515            6.12511827e-06,
516            9.51121820e-05,
517            3.74835076e-05,
518            4.97478018e-05,
519            ]
520        self['Patch2022']['Sf_pas'] = np.array([
521            51591.64007141,
522            52277.48834349,
523            22370.46712926,
524            22051.80366169,
525            23278.22158555,
526            ]) * 0.0349
527        self['Patch2022']['Sf_act'] = [84e3, 84e3, 120e3, 120e3, 120e3]
528
529        # Valves
530        self['Valve2022']['adaptation_A_open_fac'] = [
531            1.,
532            1.5,
533            1.,
534            1.,
535            1.5,
536            1.]
537        self['Valve2022']['A_open'] = [
538            0.00050078,
539            0.0007021,
540            0.00046806,
541            0.00051776,
542            0.00074755,
543            0.00049837]
544        self['Valve2022']['A_leak'] = [
545            2.64705882e-04,
546            2.64705882e-10,
547            2.64705882e-10,
548            2.64705882e-04,
549            2.64705882e-10,
550            2.64705882e-10]
551        self['Valve2022']['l'] = 0.01626978
552        self['Valve2022']['rho_b'] = 1050
553        self['Valve2022']['papillary_muscles'] = False
554        self['Valve2022']['papillary_muscles_slope'] = 100
555        self['Valve2022']['papillary_muscles_min'] = 0.1
556        self['Valve2022']['papillary_muscles_A_open_fac'] = 0.1
557        self['Valve2022']['soft_closure'] = True
558
559        self['Bag']['k'] = 10
560        self['Bag']['V_ref'] = 0.0005
561        self['Bag']['p_ref'] = 100
562
563        self.set('Model.Peri.TriSeg.Y', 0.035)
564
565        self['Timings']['law_tauAv'] = 2
566        self['Timings']['c_tauAv0'] = 0.2
567        self['Timings']['c_tauAv1'] = -0.485*60e-3
568
569        # Adaptation
570        self['Patch2022']['adapt_gamma'] = 0.75
571
572        self.run(stable=True)
573
574    def get_unittest_targets(self):
575        """Hardcoded results after initializing and running 1 beat."""
576        return {
577            'LVEDV': 141.6,
578            'LVESV':  69.4,
579            }
580
581
582class VanOsta2023_2308(VanOsta2023):
583    def set_reference(self):
584        self['Chamber2022']['buckling'] = False
585        self['Valve2022']['soft_closure'] = True
586        self['Valve2022']['papillary_muscles'] = True
587
588        self['Solver']['dt'] = 0.001
589        self['Solver']['dt_export'] = 0.002
590        self.set('Solver.order',  2)
591
592        self.set('Model.t_cycle', 0.85)
593
594        self.set('Model.PFC.fac',  0.5)
595        self.set('Model.PFC.epsilon', 0.1)
596        self.set('Model.PFC.fac_pfc', 1)
597        self.set('Model.PFC.stable_threshold', 0.001)
598
599        # Set parameters Circulation
600        self['Tube0D']['A0'] = [0.0004982 , 0.00049959, 0.00045184, 0.00051768]
601        self['Tube0D']['A_wall'] = [1.13597208e-04, 3.79106863e-05, 8.91394588e-05, 4.26698697e-05]
602        self['Tube0D']['k'] = [1.66666667, 2.33333333, 1.        , 2.33333333]
603        self['Tube0D']['l'] = [0.4, 0.4, 0.1, 0.1]
604        self['Tube0D']['p0'] = [12154.79898845,   213.26733196,  1913.05480512,   600.72294959]
605
606        self['ArtVen']['p0'] = [6345.26731406,  950.        ]
607        self['ArtVen']['q0'] = [4.5e-05, 4.5e-05]
608        self['ArtVen']['k'] = [1.  , 1.72]
609
610        # Set volume state variables
611        # not needed to reset, only if initial model crashes
612        if True:
613            self.set('Model.SyArt.V', 2e-4)
614            self.set('Model.PuArt.V', 1e-4)
615            self.set('Model.SyVen.V', 4e-4)
616            self.set('Model.PuVen.V', 2e-4)
617
618            self.set('Model.Peri.TriSeg.V', 56e-6)
619            self.set('Model.Peri.TriSeg.Y', 36e-3)
620
621            self.set('Model.Peri.TriSeg.cLv.V',  150e-6)
622            self.set('Model.Peri.TriSeg.cRv.V',  110e-6)
623            self.set('Model.Peri.La.V',  25e-6)
624            self.set('Model.Peri.Ra.V',  150e-6)
625
626        # all wall parameters
627        self['Patch2022']['l_se'] = 0.04
628        self['Patch2022']['l_s_ref'] = 2.0
629        self['Patch2022']['l_s0'] = 1.8
630        self['Patch2022']['dl_s_pas'] = 0.6
631        self['Patch2022']['k1'] = 10
632        self['Patch2022']['dt'] = 0
633        self['Patch2022']['C_rest'] = 0
634        self['Patch2022']['l_si0'] = 1.51
635        self['Patch2022']['LDAD'] = [1.057, 1.057, 0.64 , 0.64 , 0.64 ]
636        self['Patch2022']['ADO'] = [0.65, 0.65, 0.75, 0.75, 0.75]
637        self['Patch2022']['LDCC'] = [4. , 4. , 3.2, 3.2, 3.2]
638        self['Patch2022']['v_max'] = 7.
639        self['Patch2022']['v_max'][:2] = 14.
640
641        self['Patch2022']['tr'] = [0.4 , 0.4 , 0.24, 0.24, 0.24]
642        self['Patch2022']['td'] = [0.4 , 0.4 , 0.23, 0.23, 0.23]
643
644        # wall specific
645        self['Patch2022']['Sf_act'] = [80e3, 80e3, 120e3, 120e3, 120e3]
646        self['Patch2022']['Am_ref'] = [0.00458522, 0.00381765, 0.00806997, 0.00454307, 0.01176514]
647        self['Patch2022']['V_wall'] = [2.75050615e-05, 1.25959326e-05, 8.80910385e-05, 3.46281364e-05,
648               5.16086110e-05]
649        self['Patch2022']['Sf_pas'] = [ 16.51117993,  17.0639141 , 580.39235927, 569.67392472,
650               606.50514861]
651
652        # Valves
653        self['Valve2022']['adaptation_A_open_fac'] = [1.  , 1.11, 1.  , 1.  , 1.11, 1.  ]
654        self['Valve2022']['A_open'] = [0.00050042, 0.00050203, 0.00045228, 0.00051764, 0.00055321,
655               0.00049838]
656        self['Valve2022']['A_leak'] = [2.64705882e-04, 2.64705882e-10, 2.64705882e-10, 2.64705882e-04,
657               2.64705882e-10, 2.64705882e-10]
658        self['Valve2022']['l'] = 0.01626978
659        self['Valve2022']['rho_b'] = 1050
660        self['Valve2022']['papillary_muscles'] = True
661        self['Valve2022']['papillary_muscles_slope'] = 100
662        self['Valve2022']['papillary_muscles_min'] = 0.1
663        self['Valve2022']['papillary_muscles_A_open_fac'] = 0.1
664        self['Valve2022']['soft_closure'] = True
665
666        self['Bag']['k'] = 10
667        self['Bag']['V_ref'] = [0.00051968]
668        self['Bag']['p_ref'] = 100
669
670        self.set('Model.Peri.TriSeg.Y', 0.035)
671
672        self['Timings']['law_tauAv'] = 2
673        self['Timings']['c_tauAv0'] = 0.172
674        self['Timings']['c_tauAv1'] = -0.485*60e-3
675
676        # Adaptation
677        self['Patch2022']['adapt_gamma'] = 0.5
678        self.set('Solver.store_beats', 1)
679
680        self['PressureFlowControl']['stable_threshold'] = 1e-3
681
682        self['Patch2022']['SfPasMaxT'] = [6400., 6400.,  6600.,  6600.,  6600.]
683        self['Patch2022']['FacSfActT'] = [0.44, 0.44, 0.61, 0.61, 0.61]
684        self['Patch2022']['SfPasActT'] = [4800., 4800., 6600., 6600., 6600.]
685        self['Patch2022']['LsPasActT'] = [3.  , 3.  , 2.23, 2.23, 2.23]
686
687        self.run(stable=True)
688
689        return
690
691    def get_unittest_targets(self):
692        """Hardcoded results after initializing and running 1 beat."""
693        return {
694            'LVEDV': 119.7,
695            'LVESV':  47.2,
696            }