Simple One-Ventricle Open Circulation
Todo
Make One-Ventricle Circulation figure
(Source code, png, hires.png, pdf)
The full code to generate this plot is shown below.
1# -*- coding: utf-8 -*-
2import circadapt
3
4import numpy as np
5import matplotlib.pyplot as plt
6
7import time
8
9def create_model():
10 ###
11 n_beat = 100
12
13 dt = 0.001
14 solver = "forward_euler"
15
16 # open model
17 model = circadapt.CircAdapt(solver)
18
19 model.set('Solver.dt', dt)
20 model.set('Solver.dt_export', dt)
21
22 # Build the cavity C using a Chamber2022 object. This object automatically
23 # gets an Wall2022 object with no patches, so a patch must be added.
24 model.add_component('Chamber2022', 'C')
25 model.add_component('Patch2022', 'P', 'C.wC')
26
27 # add windkessel
28 model.add_component('NodePressure', 'Prox')
29 model.add_component('NodePressure', 'Dist')
30 model.add_component('Capacitor', 'Art')
31 model.add_component('Resistance', 'Resistance')
32
33 model.add_component('Diode', 'ValveIn')
34 model.add_component('Diode', 'ValveOut')
35
36 model.set_component('ValveIn.Prox', 'Prox')
37 model.set_component('ValveIn.Dist', 'C')
38
39 model.set_component('ValveOut.Prox', 'C')
40 model.set_component('ValveOut.Dist', 'Art')
41
42 model.set_component('Resistance.Prox', 'Art')
43 model.set_component('Resistance.Dist', 'Dist')
44
45 # parameters
46 # model.set('Model.Valve.R', 10)
47 model.set('Model.ValveIn.R', 10e5)
48 model.set('Model.ValveOut.R', 1e7)
49 model.set('Model.Resistance.R', 1e8)
50 model.set('Model.Art.C', 1e-8)
51
52 # set proximal pressure
53 model.set('Model.Prox.Pressure.set_time', 0)
54 model.set('Model.Prox.Pressure.set_value', 2000)
55
56
57 # set state variables
58 model.set('Model.C.V', 10e-6)
59 # model.set('Model.SyArt.V', 200e-6)
60 # model.set('Model.SyVen.V', 300e-6)
61
62 # parameterize patch
63 model['Patch2022']['dt'] = 0.1
64 model['Patch2022']['Sf_act'] = 100e3
65 model['Patch2022']['Sf_pas'] = 1e3
66 model['Patch2022']['Am_ref'] = 0.011
67 model['Patch2022']['k1'] = 10
68 model['Patch2022']['v_max'] = 7
69 model['Patch2022']['V_wall'] = 1e-04
70 model['Patch2022']['tr'] = 0.25
71 model['Patch2022']['td'] = 0.25
72 model['Patch2022']['time_act'] = 0.2
73 model.set('Model.C.wC.P.l_si', 2.0)
74 model.set('Model.Art.V', 10e-4)
75 model.set('Model.Art.V0', 1e-4)
76 # model.set('Model.Art.C', 2e-8)
77 # model.set('Model.Art.V', 200)
78
79
80 # model['Tube0D']['p0'][0] = 900
81 # model['Tube0D']['A0'][0] = 6e-4
82
83 # Run beats
84 t0 = time.time()
85 model.run(n_beat)
86 print(time.time()-t0)
87
88 return model
89
90# %% Run model and plot
91if __name__ == '__main__':
92 model = create_model()
93 plt.figure(1, clear=True, figsize=(12, 4))
94
95 t = model.get('Solver.t') * 1e3
96
97 m=1
98 n=3
99
100 ax = plt.subplot(m,n,1)
101 ax.plot(t, model.get('Model.C.p')/133, label='Chamber', color='k')
102 ax.plot(t, model.get('Model.Prox.p')/133, label='Prox', color='r')
103 ax.plot(t, model.get('Model.Art.p')/133, label='Art', color='r')
104 ax.plot(t, model.get('Model.Dist.p')/133, label='Dist', color='b')
105 plt.legend()
106 ax.set_ylabel('Pressure [mmHg]')
107 ax.set_xlabel('Time [ms]')
108
109 ax = plt.subplot(m,n,2)
110 ax.plot(t, model.get('Model.C.V')*1e6, label='Chamber', color='k')
111 ax.set_ylabel('Volume [mL]')
112 ax.set_xlabel('Time [ms]')
113 plt.legend()
114
115 ax = plt.subplot(m,n,3)
116 ax.plot(t, model.get('Model.ValveIn.q')*1e3, label='ValveIn')
117 ax.plot(t, model.get('Model.ValveOut.q')*1e3, label='ValveOut')
118 ax.plot(t, model.get('Model.Resistance.q')*1e3, label='Resistance')
119 plt.legend()
120
121 ax.axhline(0, color='k', linestyle='--')
122 ax.set_ylabel('Flow [mL/ms]')
123 ax.set_xlabel('Time [ms]')
124 plt.legend()
125
126 plt.suptitle('Simple one-ventricle model setup', fontsize=16, weight='bold')
127
128
129 # ax = plt.subplot(m,n,4)
130 # ax.axhline(0, color='k', linestyle='--', lw=1)
131 # ax.set_ylabel('C')
132 # plt.legend()
133
134 plt.tight_layout()
135 plt.draw()
136 plt.show()
137