Preload Afterload experiment
Tutorial on how to implement the preload-afterload experiment. This tutorial discusses how to build the model, how to switch between different patch types, and how to extract and plot information from the model.
Todo
Make Preload Afterload methods figure
As all projects, we start with importing packages.
python
1from circadapt import CircAdapt
2import numpy as np
3import matplotlib.pyplot as plt
4
5import time
6
7
We then create an empty model. Then, we add the PreAfterloadExperiment wall object. This object is empty by default, but it needs at least 1 patch object to work. In this tutorial, we use the Patch2022 object.
python
1model = CircAdapt(
2 'Custom',
3 'SolverFE',
4 )
5model.add_component('PreAfterloadExperiment', 'PAE')
6model.add_component('Patch2022', 'P', 'PAE')
7
The model has to be parametized, as default parameter values do not make sense for this setup.
python
1model.set('Model.tCycle', 1e-0)
2model.set('Solver.dT', 1e-3)
3model.set('Solver.dTexport', 1e-3)
4
5model['PreAfterloadExperiment']['AmRefAfterload'] = 0.009
6model['PreAfterloadExperiment']['Tafterload'] = 200
7model['PreAfterloadExperiment']['nIter'] = 5
8model['Patch2022']['AmRef'] = 0.0076
9model['Patch2022']['VWall'] = 92.43*1e-6
10model['Patch2022']['ADO'] = 3.0
11model['Patch2022']['TR'] = 0.5
12model['Patch2022']['TD'] = 0.5
13model['Patch2022']['k1'] = 10.
14model['Patch2022']['vMax'] = 7.0
15
16model.set('Model.PAE.P.Lsi', -0.04 + 2 * np.sqrt(
17 model['PreAfterloadExperiment']['AmRefAfterload']/model['Patch2022']['AmRef']))
18
Now the model setup is finished and we can run the simulation. For this setup, only 1 run is sufficient.
python
1t0 = time.time()
2model.run(1)
3t1 = time.time()
4print(t1-t0)
5
The simulation will result in the stresses and tensions shown in the plot below.
(Source code, png, hires.png, pdf)
The full code is shown below.
python
1from circadapt import CircAdapt
2import numpy as np
3import matplotlib.pyplot as plt
4
5import time
6
7
8# %% Create custom model
9model = CircAdapt(
10 'Custom',
11 'SolverFE',
12 )
13model.add_component('PreAfterloadExperiment', 'PAE')
14model.add_component('Patch2022', 'P', 'PAE')
15
16# %% Set model parameters
17model.set('Model.tCycle', 1e-0)
18model.set('Solver.dT', 1e-3)
19model.set('Solver.dTexport', 1e-3)
20
21model['PreAfterloadExperiment']['AmRefAfterload'] = 0.009
22model['PreAfterloadExperiment']['Tafterload'] = 200
23model['PreAfterloadExperiment']['nIter'] = 5
24model['Patch2022']['AmRef'] = 0.0076
25model['Patch2022']['VWall'] = 92.43*1e-6
26model['Patch2022']['ADO'] = 3.0
27model['Patch2022']['TR'] = 0.5
28model['Patch2022']['TD'] = 0.5
29model['Patch2022']['k1'] = 10.
30model['Patch2022']['vMax'] = 7.0
31
32model.set('Model.PAE.P.Lsi', -0.04 + 2 * np.sqrt(
33 model['PreAfterloadExperiment']['AmRefAfterload']/model['Patch2022']['AmRef']))
34
35# %% Run model
36t0 = time.time()
37model.run(1)
38t1 = time.time()
39print(t1-t0)
40
41# %% Plot model
42fig = plt.figure(1, clear=True, figsize=(12, 8))
43m = 2
44n = 3
45
46t = model.get('Solver.Time') * 1e3
47
48ax1 = fig.add_subplot(m, n, 1)
49ax1.plot(t, model.get('Model.PAE.P.Ls'))
50ax1.plot(t, model.get('Model.PAE.P.Lsi'))
51ax1.legend(['Ls', 'Lsi'])
52ax1.set_ylabel('Sarcomere length [$\mu m$]')
53
54ax4 = fig.add_subplot(m, n, 4)
55ax4.plot(t, model.get('Model.PAE.P.Am')*1e4, label='Patch Am')
56ax4.plot(t, model.get('Model.PAE.P.Am0')*1e4, label='Patch Am0')
57ax4.set_ylabel('Segment area [$cm^2$]')
58ax4.legend()
59
60ax3 = fig.add_subplot(m, n, 3)
61ax3.plot(t, model.get('Model.PAE.P.Sf')*1e-3, label='Sf')
62ax3.plot(t, model.get('Model.PAE.P.SfEcm')*1e-3, label='SfEcm')
63ax3.set_ylabel('Stress [kPa]')
64ax3.legend()
65
66ax2 = fig.add_subplot(m, n, 2)
67ax2.plot(t, model.get('Model.PAE.P.C'), label='C')
68ax2.set_ylabel('C-curve [-]')
69
70ax5 = fig.add_subplot(m, n, 5)
71ax5.plot(t, model.get('Model.PAE.DADT')*1e4, label='Tension Wall')
72ax5.set_ylabel('DADT [$cm^2/N$]')
73
74ax6 = fig.add_subplot(m, n, 6)
75ax6.plot(t, model.get('Model.PAE.T'), label='Tension Wall')
76ax6.set_ylabel('Tension [N]')
77
78# Plot design
79for ax in [ax1, ax2, ax3, ax4, ax5, ax6]:
80 ax.spines[['right', 'top']].set_visible(False)
81ax5.set_xlabel('Time [ms]')
82
83plt.tight_layout()