GitHub

Mechanics Tutorial 3: Coupling with Lumped Blood Circuits

Pressure Volume Loop

This tutorial shows how to couple 3d chamber models with 0d fluid models.

Introduction

In this tutorial we will reproduce a simplified version of the model presented by Regazzoni et al. [3].

Warning

The API for 3D-0D coupling is work in progress and is hence subject to potential breaking changes.

Commented Program

We start by loading Thunderbolt and LinearSolve to use a custom direct solver of our choice.

using Thunderbolt, LinearSolve

Finally, we try to approach a valid initial state by solving a simpler model first.

using OrdinaryDiffEqTsit5, OrdinaryDiffEqOperatorSplitting

fluid_model_init = RSAFDQ2022LumpedCicuitModel()
u0 = zeros(Thunderbolt.num_states(fluid_model_init))
Thunderbolt.default_initial_condition!(u0, fluid_model_init)
prob = ODEProblem((du, u, p, t) -> Thunderbolt.lumped_driver!(du, u, t, [], p), u0, (0.0, 100*fluid_model_init.THB), fluid_model_init)
sol = solve(prob, Tsit5())
retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 11931-element Vector{Float64}:
     0.0
     0.04838838098003728
     0.2792640643984848
     0.7267523496016612
     1.3438604381655068
     2.144258889210636
     3.155553117293995
     4.398246405764439
     5.905268886754168
     7.716402419487973
     ⋮
 79949.53923019694
 79956.29777079161
 79963.05665516232
 79969.81596119265
 79976.57556892763
 79983.3352996462
 79990.09501501129
 79996.8546516533
 80000.0
u: 11931-element Vector{Vector{Float64}}:
 [65.0, 120.0, 65.0, 145.0, 10.66, 4.0, 4.67, 3.2, 0.0, 0.0, 0.0, 0.0]
 [64.99775599207983, 120.00227770768569, 64.99411468494874, 145.00593718603966, 10.659998698778578, 3.9999999665140544, 4.669999657362968, 3.199999933568232, 0.0004815730430550299, 0.002134255662419863, 0.001058704949551402, 0.0013873724257603324]
 [64.98801864146336, 120.01307519006716, 64.96751009486455, 145.03415042166324, 10.65995732473332, 3.9999989362827226, 4.66998887341767, 3.1999978409188183, 0.0027286563340648087, 0.011610931373988172, 0.005887067909681074, 0.007715375293357945]
 [64.97343673901082, 120.03362955227148, 64.92214626693816, 145.08829571065598, 10.659717831193655, 3.9999934149906573, 4.6699281242379485, 3.1999860448630866, 0.006854505222657548, 0.027039465541638375, 0.014273965931167643, 0.018710194902465015]
 [64.96169085679992, 120.06099336358118, 64.87058951967002, 145.161795531686, 10.659065710108699, 3.9999800212165577, 4.669769295052095, 3.19995517731026, 0.012080360027470989, 0.0432139565112982, 0.024009535416236124, 0.03147877059461845]
 [64.9590422884215, 120.0943385463058, 64.81815518830922, 145.25505063752837, 10.657716529478249, 3.999956125179102, 4.669457198522932, 3.1998944573339676, 0.018132298859120682, 0.057781487036428374, 0.03404823736490992, 0.0446530382982922]
 [64.97315242052815, 120.13201303356213, 64.76867906196281, 145.36931403553868, 10.655297467403598, 3.99992028271473, 4.668930854853903, 3.199791923888736, 0.024747263721687698, 0.06934017382066814, 0.043508831911020375, 0.05707863744246779]
 [65.01361095019732, 120.16952228339345, 64.72592307611659, 145.50371995376398, 10.651396546185397, 3.999873170375802, 4.668138641692197, 3.1996373849033297, 0.03153312463200297, 0.07731247060104394, 0.051579011162919436, 0.06768823773954884]
 [65.09351105121917, 120.1983518552788, 64.6935907758836, 145.65663985325492, 10.645541833828794, 3.9998165986079437, 4.667034206177342, 3.1994216580901362, 0.038132755971079714, 0.08210581363088908, 0.05784520570133698, 0.07593295740565024]
 [65.2293722454306, 120.20576277678236, 64.67799794660367, 145.8235332540356, 10.637227798674608, 3.9997528795818074, 4.665579203954888, 3.1991372068563617, 0.04420483062644496, 0.084592086707584, 0.06223503408143795, 0.0817052270314148]
 ⋮
 [63.19267337105863, 114.00416580856115, 60.497952220122066, 146.6449554353661, 11.69189003308238, 3.9910356533134648, 4.88046481280082, 3.169539343878219, 0.07267704475187975, 0.08702602805204065, 0.07911179376397252, 0.08824125124438174]
 [63.1841831374728, 114.60940124917967, 60.680747887776974, 147.0507811537661, 11.637509612609808, 3.9909530253037753, 4.873350062357445, 3.1690136319353135, 0.07216457877707018, 0.0869993148634757, 0.07880657057600357, 0.08833815794614759]
 [63.229074017074254, 115.1617785321413, 60.876243914366285, 147.44376596512296, 11.583509893988257, 3.9908676591723036, 4.866262349430057, 3.168466663657466, 0.07165573158473659, 0.08697571033200799, 0.07850364517224674, 0.0883888939025456]
 [63.31839645943162, 115.66997727483837, 61.0846248665232, 147.82375306078924, 11.529887592387322, 3.9907795564250454, 4.859201394843491, 3.1679008627909524, 0.0711504721612215, 0.08695515109188282, 0.07820287909975845, 0.08840350887344126]
 [63.44499195480075, 116.14092948524322, 61.30605031708183, 148.19059302553396, 11.476641011705432, 3.9906887213184654, 4.852167137023105, 3.16731815518271, 0.07064878446129669, 0.08693759749933908, 0.07790416966598689, 0.08839006292355864]
 [63.60313545620179, 116.58017737863477, 61.54066041767384, 148.54415045277733, 11.423768905459816, 3.990595159247043, 4.845159581986698, 3.1667200605463917, 0.0701506566722406, 0.08692302112177575, 0.07760743822453335, 0.08835501976231346]
 [63.7882520892354, 116.99215576791856, 61.78858112374743, 148.88430740071777, 11.371269693661098, 3.9904988753829604, 4.838178698274948, 3.166107765629687, 0.06965607383539181, 0.08691139577280431, 0.07731262150798661, 0.08830356026249267]
 [63.99669032244723, 117.38041416755063, 62.04992764998896, 149.21096297654444, 11.319141208644952, 3.9903998740495217, 4.831224380872002, 3.1654821861418836, 0.06916501545320265, 0.08690269378070858, 0.07701966675734155, 0.0882398340563524]
 [64.10079198906878, 117.55380932929918, 62.176016658822775, 149.35833811202073, 11.295010573762173, 3.990352984138106, 4.827997427577099, 3.165186764819986, 0.06893771251040942, 0.08696357980372327, 0.07688397157456074, 0.08820688785869524]

plot(sol, idxs=[1,2,3,4], tspan=(99fluidmodelinit.THB, 100fluidmodelinit.THB))

# Precomputed initial guess
u₀fluid = sol.u[end]
@info "Total blood volume: $(sum(u₀fluid[1:4])) + $(fluid_model_init.Csysₐᵣ*u₀fluid[5]) + $(fluid_model_init.Csysᵥₑₙ*u₀fluid[6]) + $(fluid_model_init.Cpulₐᵣ*u₀fluid[7]) + $(fluid_model_init.Cpulᵥₑₙ*u₀fluid[8])"
[ Info: Total blood volume: 393.1889560892115 + 101.66345564540934 + 4788.81739818597 + 362.1295874619882 + 379.8536497430151

We now generate the mechanical subproblem as in the first tutorial

scaling_factor = 3.7;
Warning

Tuning parameter until all bugs are fixed in this tutorial :)

mesh = generate_ideal_lv_mesh(8,2,5;
    inner_radius = scaling_factor*0.7,
    outer_radius = scaling_factor*1.0,
    longitudinal_upper = 0.4,
    apex_inner = scaling_factor* 1.3,
    apex_outer = scaling_factor*1.5
)
mesh = Thunderbolt.hexahedralize(mesh)
Thunderbolt.SimpleMesh{3, Hexahedron, Float64} with 736 Hexahedron cells and 965 nodes
  Surface subdomains:
    Epicardium 184 Hexahedron cells
    Base 64 Hexahedron cells
    Endocardium 184 Hexahedron cells
Todo

The 3D0D coupling does not yet support multiple subdomains.

coordinate_system = compute_lv_coordinate_system(mesh)
microstructure    = create_microstructure_model(
    coordinate_system,
    LagrangeCollection{1}()^3,
    ODB25LTMicrostructureParameters(),
);
passive_material_model = Guccione1991PassiveModel()
active_material_model  = Guccione1993ActiveModel()
function calcium_profile_function(x::LVCoordinate,t_global)
    linear_interpolation(t,y1,y2,t1,t2) = y1 + (t-t1) * (y2-y1)/(t2-t1)
    ca_peak(x)                          = 1.0
    t = t_global % 800.0
    if 0 ≤ t ≤ 120.0
        return linear_interpolation(t,        0.0, ca_peak(x),   0.0, 120.0)
    elseif t ≤ 272.0
        return linear_interpolation(t, ca_peak(x),        0.0, 120.0, 272.0)
    else
        return 0.0
    end
end
calcium_field = AnalyticalCoefficient(
    calcium_profile_function,
    coordinate_system,
)
sarcomere_model = CaDrivenInternalSarcomereModel(ConstantStretchModel(), calcium_field)
active_stress_model = ActiveStressModel(
    passive_material_model,
    active_material_model,
    sarcomere_model,
    microstructure,
)
weak_boundary_conditions = (RobinBC(1.0, "Epicardium"),NormalSpringBC(100.0, "Base"))
solid_model = QuasiStaticModel(:displacement, active_stress_model, weak_boundary_conditions);

The solid model is now couple with the circuit model by adding a Lagrange multipliers constraining the 3D chamber volume to match the chamber volume in the 0D model.

fluid_model = RSAFDQ2022LumpedCicuitModel(; lv_pressure_given = false)
coupler = LumpedFluidSolidCoupler(
    [
        ChamberVolumeCoupling(
            "Endocardium",
            RSAFDQ2022SurrogateVolume(),
            :Vₗᵥ,
            :pₗᵥ,
        )
    ],
    :displacement,
)
coupled_model = RSAFDQ2022Model(solid_model, fluid_model, coupler);
Todo

Once we figure out a nicer way to do this we should add more detailed docs here.

Now we semidiscretize the model spatially as usual with finite elements and annotate the model with a stable split.

spatial_discretization_method = FiniteElementDiscretization(
    Dict(:displacement => LagrangeCollection{1}()^3),
    [
        Dirichlet(:displacement, getfacetset(mesh, "Base"), (x,t) -> [0.0], [3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor1"), (x,t) -> (0.0, 0.0, 0.0), [1,2,3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor2"), (x,t) -> (0.0, 0.0), [2,3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor3"), (x,t) -> (0.0,), [3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor4"), (x,t) -> (0.0,), [3])
    ],
)
splitform = semidiscretize(
    RSAFDQ2022Split(coupled_model),
    spatial_discretization_method,
    mesh,
)

dt₀ = 1.0
dtvis = 10.0
tspan = (0.0, 3*800.0)
(0.0, 2400.0)

This speeds up the CI # hide

(0.0, 10.0)

The remaining code is very similar to how we use SciML solvers.

chamber_solver = HomotopyPathSolver(
    NewtonRaphsonSolver(;
        max_iter=10,
        tol=1e-2,
        inner_solver=SchurComplementLinearSolver(
            LinearSolve.UMFPACKFactorization()
        )
    )
)
blood_circuit_solver = Tsit5()
timestepper = LieTrotterGodunov((chamber_solver, blood_circuit_solver))

u₀ = zeros(solution_size(splitform))
u₀solid_view = @view  u₀[OS.get_solution_indices(splitform, 1)]
u₀fluid_view = @view  u₀[OS.get_solution_indices(splitform, 2)]
u₀fluid_view .= u₀fluid

problem = OperatorSplittingProblem(splitform, u₀, tspan)
integrator = init(problem, timestepper, dt=dt₀, verbose=true; dtmax=10.0);

# f2 = Figure()
# axs = [
#     Axis(f2[1, 1], title="LV"),
#     Axis(f2[1, 2], title="RV"),
#     Axis(f2[2, 1], title="LA"),
#     Axis(f2[2, 2], title="RA")
# ]

# vlv = Observable(Float64[])
# plv = Observable(Float64[])

# vrv = Observable(Float64[])
# prv = Observable(Float64[])

# vla = Observable(Float64[])
# pla = Observable(Float64[])

# vra = Observable(Float64[])
# pra = Observable(Float64[])

# lines!(axs[1], vlv, plv)
# lines!(axs[2], vrv, prv)
# lines!(axs[3], vla, pla)
# lines!(axs[4], vra, pra)
# for i in 1:4
#     xlims!(axs[1], 0.0, 180.0)
#     ylims!(axs[1], 0.0, 180.0)
# end
# display(f2)
t: 0.0
u: 2908-element Vector{Float64}:
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   ⋮
 149.35833811202073
  11.295010573762173
   3.990352984138106
   4.827997427577099
   3.165186764819986
   0.06893771251040942
   0.08696357980372327
   0.07688397157456074
   0.08820688785869524
Todo

recover online visualization of the pressure volume loop

Todo

The post-processing API is not yet finished. Please revisit the tutorial later to see how to post-process the simulation online. Right now the solution is just exported into VTK, such that users can visualize the solution in e.g. ParaView.

Now we can finally solve the coupled problem in time.

io = ParaViewWriter("CM03_3d0d-coupling");
for (u, t) in TimeChoiceIterator(integrator, tspan[1]:dtvis:tspan[2])
    chamber_function = OS.get_operator(splitform, 1)
    (; dh) = chamber_function.structural_function
    store_timestep!(io, t, dh.grid)
    usolid_view = @view u[OS.get_solution_indices(splitform, 1)]
    Thunderbolt.store_timestep_field!(io, t, dh, usolid_view, :displacement)
    Thunderbolt.finalize_timestep!(io, t)

    # if t > 0.0
    #     lv = chamber_function.tying_info.chambers[1]
    #     append!(vlv.val, lv.V⁰ᴰval)
    #     append!(plv.val, u[lv.pressure_dof_index_global])
    #     notify(vlv)
    #     notify(plv)
    # end
    # TODO plot other chambers
end
[ Info: Volume 3D 112.23695369430204
[ Info: Chamber 1 p=0.0, V0=117.55380932929918
[ Info: Volume 3D 117.62016991008453
[ Info: Chamber 1 p=0.3233629322229942, V0=117.55380932929918
[ Info: Volume 3D 117.55381416933325
[ Info: Chamber 1 p=0.3177275538568636, V0=117.55380932929918
[ Info: Volume 3D 117.55380907641506
[ Info: Chamber 1 p=0.31774697598283924, V0=117.55380932929918
[ Info: Volume 3D 117.55380907641506
[ Info: Chamber 1 p=0.31774697598283924, V0=118.48189820240242
[ Info: Volume 3D 118.48735051246099
[ Info: Chamber 1 p=0.5585722931067562, V0=118.48189820240242
[ Info: Volume 3D 118.48191739099943
[ Info: Chamber 1 p=0.558272116204863, V0=118.48189820240242
[ Info: Volume 3D 118.48191739099943
[ Info: Chamber 1 p=0.558272116204863, V0=119.15406122935076
[ Info: Volume 3D 119.15813146506197
[ Info: Chamber 1 p=0.7932402092818922, V0=119.15406122935076
[ Info: Volume 3D 119.15407803529165
[ Info: Chamber 1 p=0.7930988857578978, V0=119.15406122935076
[ Info: Volume 3D 119.15407803529165
[ Info: Chamber 1 p=0.7930988857578978, V0=119.58123372945208
[ Info: Volume 3D 119.58401838922957
[ Info: Chamber 1 p=1.0219980745501422, V0=119.58123372945208
[ Info: Volume 3D 119.58124583507879
[ Info: Chamber 1 p=1.022048166848415, V0=119.58123372945208
[ Info: Volume 3D 119.58124583507879
[ Info: Chamber 1 p=1.022048166848415, V0=119.77429689099341
[ Info: Volume 3D 119.77623460309154
[ Info: Chamber 1 p=1.2449025622928083, V0=119.77429689099341
[ Info: Volume 3D 119.77430546718797
[ Info: Chamber 1 p=1.2451071184912486, V0=119.77429689099341
[ Info: Volume 3D 119.77430546718797
[ Info: Chamber 1 p=1.2451071184912486, V0=119.77429788957461
[ Info: Volume 3D 119.77578062886978
[ Info: Chamber 1 p=1.462562397912275, V0=119.77429788957461
[ Info: Volume 3D 119.77430409906465
[ Info: Chamber 1 p=1.4628717432198004, V0=119.77429788957461
[ Info: Volume 3D 119.77430409906465
[ Info: Chamber 1 p=1.4628717432198004, V0=119.7742988439663
[ Info: Volume 3D 119.77555784761024
[ Info: Chamber 1 p=1.6789513610207192, V0=119.7742988439663
[ Info: Volume 3D 119.77430338386469
[ Info: Chamber 1 p=1.6792588369140777, V0=119.7742988439663
[ Info: Volume 3D 119.77430338386469
[ Info: Chamber 1 p=1.6792588369140777, V0=119.77429975444441
[ Info: Volume 3D 119.7753863829272
[ Info: Chamber 1 p=1.8940857797181911, V0=119.77429975444441
[ Info: Volume 3D 119.77430319612104
[ Info: Chamber 1 p=1.8943885853517342, V0=119.77429975444441
[ Info: Volume 3D 119.77430319612104
[ Info: Chamber 1 p=1.8943885853517342, V0=119.77430062126072
[ Info: Volume 3D 119.775248997199
[ Info: Chamber 1 p=2.1080639354756023, V0=119.77430062126072
[ Info: Volume 3D 119.77430330819017
[ Info: Chamber 1 p=2.1083610623028717, V0=119.77430062126072
[ Info: Volume 3D 119.77430330819017
[ Info: Chamber 1 p=2.1083610623028717, V0=119.77430144464691
[ Info: Volume 3D 119.77513576558773
[ Info: Chamber 1 p=2.3209701087149703, V0=119.77430144464691
[ Info: Volume 3D 119.77430359255796
[ Info: Chamber 1 p=2.3212613596851703, V0=119.77430144464691
Tip

If you want to see more details of the solution process launch Julia with Thunderbolt as debug module:

JULIA_DEBUG=Thunderbolt julia --project --threads=auto my_simulation_runner.jl

References

[3]
F. Regazzoni, M. Salvador, P. C. Africa, M. Fedele, L. Dedè and A. Quarteroni. A cardiac electromechanical model coupled with a lumped-parameter model for closed-loop blood circulation. Journal of Computational Physics 457, 111083 (2022).

Plain program

Here follows a version of the program without any comments. The file is also available here: cm03_3d0d-coupling.jl.

using Thunderbolt, LinearSolve

using OrdinaryDiffEqTsit5, OrdinaryDiffEqOperatorSplitting

fluid_model_init = RSAFDQ2022LumpedCicuitModel()
u0 = zeros(Thunderbolt.num_states(fluid_model_init))
Thunderbolt.default_initial_condition!(u0, fluid_model_init)
prob = ODEProblem((du, u, p, t) -> Thunderbolt.lumped_driver!(du, u, t, [], p), u0, (0.0, 100*fluid_model_init.THB), fluid_model_init)
sol = solve(prob, Tsit5())

# Precomputed initial guess
u₀fluid = sol.u[end]
@info "Total blood volume: $(sum(u₀fluid[1:4])) + $(fluid_model_init.Csysₐᵣ*u₀fluid[5]) + $(fluid_model_init.Csysᵥₑₙ*u₀fluid[6]) + $(fluid_model_init.Cpulₐᵣ*u₀fluid[7]) + $(fluid_model_init.Cpulᵥₑₙ*u₀fluid[8])"

scaling_factor = 3.7;

mesh = generate_ideal_lv_mesh(8,2,5;
    inner_radius = scaling_factor*0.7,
    outer_radius = scaling_factor*1.0,
    longitudinal_upper = 0.4,
    apex_inner = scaling_factor* 1.3,
    apex_outer = scaling_factor*1.5
)
mesh = Thunderbolt.hexahedralize(mesh)

coordinate_system = compute_lv_coordinate_system(mesh)
microstructure    = create_microstructure_model(
    coordinate_system,
    LagrangeCollection{1}()^3,
    ODB25LTMicrostructureParameters(),
);
passive_material_model = Guccione1991PassiveModel()
active_material_model  = Guccione1993ActiveModel()
function calcium_profile_function(x::LVCoordinate,t_global)
    linear_interpolation(t,y1,y2,t1,t2) = y1 + (t-t1) * (y2-y1)/(t2-t1)
    ca_peak(x)                          = 1.0
    t = t_global % 800.0
    if 0 ≤ t ≤ 120.0
        return linear_interpolation(t,        0.0, ca_peak(x),   0.0, 120.0)
    elseif t ≤ 272.0
        return linear_interpolation(t, ca_peak(x),        0.0, 120.0, 272.0)
    else
        return 0.0
    end
end
calcium_field = AnalyticalCoefficient(
    calcium_profile_function,
    coordinate_system,
)
sarcomere_model = CaDrivenInternalSarcomereModel(ConstantStretchModel(), calcium_field)
active_stress_model = ActiveStressModel(
    passive_material_model,
    active_material_model,
    sarcomere_model,
    microstructure,
)
weak_boundary_conditions = (RobinBC(1.0, "Epicardium"),NormalSpringBC(100.0, "Base"))
solid_model = QuasiStaticModel(:displacement, active_stress_model, weak_boundary_conditions);

fluid_model = RSAFDQ2022LumpedCicuitModel(; lv_pressure_given = false)
coupler = LumpedFluidSolidCoupler(
    [
        ChamberVolumeCoupling(
            "Endocardium",
            RSAFDQ2022SurrogateVolume(),
            :Vₗᵥ,
            :pₗᵥ,
        )
    ],
    :displacement,
)
coupled_model = RSAFDQ2022Model(solid_model, fluid_model, coupler);

spatial_discretization_method = FiniteElementDiscretization(
    Dict(:displacement => LagrangeCollection{1}()^3),
    [
        Dirichlet(:displacement, getfacetset(mesh, "Base"), (x,t) -> [0.0], [3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor1"), (x,t) -> (0.0, 0.0, 0.0), [1,2,3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor2"), (x,t) -> (0.0, 0.0), [2,3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor3"), (x,t) -> (0.0,), [3]),
        Dirichlet(:displacement, getnodeset(mesh, "MyocardialAnchor4"), (x,t) -> (0.0,), [3])
    ],
)
splitform = semidiscretize(
    RSAFDQ2022Split(coupled_model),
    spatial_discretization_method,
    mesh,
)

dt₀ = 1.0
dtvis = 10.0
tspan = (0.0, 3*800.0)

tspan = (0.0, 10.0)    # hide

chamber_solver = HomotopyPathSolver(
    NewtonRaphsonSolver(;
        max_iter=10,
        tol=1e-2,
        inner_solver=SchurComplementLinearSolver(
            LinearSolve.UMFPACKFactorization()
        )
    )
)
blood_circuit_solver = Tsit5()
timestepper = LieTrotterGodunov((chamber_solver, blood_circuit_solver))

u₀ = zeros(solution_size(splitform))
u₀solid_view = @view  u₀[OS.get_solution_indices(splitform, 1)]
u₀fluid_view = @view  u₀[OS.get_solution_indices(splitform, 2)]
u₀fluid_view .= u₀fluid

problem = OperatorSplittingProblem(splitform, u₀, tspan)
integrator = init(problem, timestepper, dt=dt₀, verbose=true; dtmax=10.0);

# f2 = Figure()
# axs = [
#     Axis(f2[1, 1], title="LV"),
#     Axis(f2[1, 2], title="RV"),
#     Axis(f2[2, 1], title="LA"),
#     Axis(f2[2, 2], title="RA")
# ]

# vlv = Observable(Float64[])
# plv = Observable(Float64[])

# vrv = Observable(Float64[])
# prv = Observable(Float64[])

# vla = Observable(Float64[])
# pla = Observable(Float64[])

# vra = Observable(Float64[])
# pra = Observable(Float64[])

# lines!(axs[1], vlv, plv)
# lines!(axs[2], vrv, prv)
# lines!(axs[3], vla, pla)
# lines!(axs[4], vra, pra)
# for i in 1:4
#     xlims!(axs[1], 0.0, 180.0)
#     ylims!(axs[1], 0.0, 180.0)
# end
# display(f2)

io = ParaViewWriter("CM03_3d0d-coupling");
for (u, t) in TimeChoiceIterator(integrator, tspan[1]:dtvis:tspan[2])
    chamber_function = OS.get_operator(splitform, 1)
    (; dh) = chamber_function.structural_function
    store_timestep!(io, t, dh.grid)
    usolid_view = @view u[OS.get_solution_indices(splitform, 1)]
    Thunderbolt.store_timestep_field!(io, t, dh, usolid_view, :displacement)
    Thunderbolt.finalize_timestep!(io, t)

    # if t > 0.0
    #     lv = chamber_function.tying_info.chambers[1]
    #     append!(vlv.val, lv.V⁰ᴰval)
    #     append!(plv.val, u[lv.pressure_dof_index_global])
    #     notify(vlv)
    #     notify(plv)
    # end
    # TODO plot other chambers
end

This page was generated using Literate.jl.