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