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