fix minilimo naming and add dynamic coupling case

Author: marinlauber

minilimo/README.md

@@ -1,11 +1,32 @@
 ## miniLIMO setups
 
-1. `limo_dynamic_coupled.jl`: works, dynamic morphing (damped) and then static Lie-Trotter coupling (weak) with the WK3 model.Pact max reaches 400 mmHg.
-2. `limo_dynamic_kostas_single.jl`: static inflation for Kostas' project. Runs well via command line Pact arguments
-3. `limo_dynamic_coupled_transient.jl`: dynamic morphing (damped) and then dynamic Lie-Trotter coupling (weak) with the WK3 model.
-4. `limo_dynamic_coupled_strong.jl`:  dynamic morphing (damped) and then dynamic, stronlgy coupled Lie-Trotter coupling with the WK3 model.
-5. `limo_dynamic_full.jl` same as the `1`, but with the full geometry, no symmetry plane.
-6. `utils`: store some repeated functions that are used across these files
+### 3D–0D Windkessel-coupled beat
+
+The four coupled scripts are named `limo_coupled_<static|dynamic>_<weak|strong>.jl`, spanning a
+2×2 design space — **structure inertia** (static equilibrium vs. dynamic HHT-α with `M·ä`) ×
+**coupling strength** (weak vs. strong):
+
+|                    | **weak** (Lie–Trotter split, black-box ODE integrator, Plv/Vlv lagged one substep) | **strong** (monolithic Newton on `(u, Plv, Pa, Pv)`, embedded implicit-Euler 0D) |
+|--------------------|---|---|
+| **static** structure  | `limo_coupled_static_weak.jl`  | `limo_coupled_static_strong.jl` |
+| **dynamic** structure | `limo_coupled_dynamic_weak.jl` | `limo_coupled_dynamic_strong.jl` |
+
+All four share the Phase-1 dynamic HHT-α morph (damped) and reload the morphed state from
+`limo_dynamic_coupled_u0.jld2`; they differ only in the Phase-2 coupled beat.
+
+1. `limo_coupled_static_weak.jl`: works, static weakly-coupled (Lie–Trotter) beat. Pact max reaches 400 mmHg. Carries `S_min`/`δp` diagnostics (Schur-complement collapse under refinement).
+2. `limo_coupled_dynamic_weak.jl`: dynamic weakly-coupled (Lie–Trotter) beat; inertia + damping retained in Phase 2.
+3. `limo_coupled_dynamic_strong.jl`: dynamic strongly (monolithically) coupled beat; HHT-α structure + implicit-Euler 0D as one Newton per step.
+4. `limo_coupled_static_strong.jl`: static counterpart of `3` — quasi-static equilibrium + implicit-Euler 0D as one Newton (drops `M·ä`/damping, `v2 = K_eff⁻¹F_plv` with no `1−α` factor).
+
+> Note: Godunov splitting ≡ Lie–Trotter. The genuine weak↔strong middle ground would be a
+> sub-iterated partitioned scheme (fixed-point Plv↔Vlv, optionally under-relaxed/Aitken) — not yet implemented.
+
+### Other setups
+
+5. `limo_dynamic_kostas_single.jl`: static inflation for Kostas' project. Runs well via command line Pact arguments.
+6. `limo_dynamic_full.jl`: same as `limo_coupled_static_weak`, but with the full geometry, no symmetry plane.
+7. `util.jl`: shared functions (mesh, edge-morph IC, RM assembly helpers) used across these files.
 
 ## old files
 

minilimo/limo_dynamic_coupled_strong.jl minilimo/limo_coupled_dynamic_strong.jlminilimo/limo_dynamic_coupled_strong.jl renamed to minilimo/limo_coupled_dynamic_strong.jl

@@ -3,7 +3,7 @@ include(joinpath(@__DIR__, "util.jl"))
 
 # Strongly (monolithically) coupled dynamic miniLIMO: HHT-α structure + 0D Windkessel
 # solved together in ONE Newton system per time step.  Counterpart of the weakly-coupled
-# `limo_dynamic_coupled_transient.jl` (Lie–Trotter split with a black-box ODE integrator).
+# `limo_coupled_dynamic_weak.jl` (Lie–Trotter split with a black-box ODE integrator).
 #
 #   PHASE 1 — dynamic morph  (t ∈ [0, T_sim], HHT-α)   [identical to the transient file]
 #     Morph the edge onto the elliptic arc + fill Plv → morphed state u = un, v, a, Plv=p_max.
@@ -194,12 +194,12 @@ u = zeros(N_dof); apply!(u, ch)
 v = zeros(N_dof)
 a = zeros(N_dof)
 
-pvd = paraview_collection("minilimo-dynamic-coupled-strong")
+pvd = paraview_collection("minilimo-coupled-dynamic-strong")
 vtk_step = Ref(0)
 resu = zeros(3, getnnodes(dh.grid))
 resθ = zeros(2, getnnodes(dh.grid))
 d, G3 = director_field(dh, scv, u)
-VTKGridFile("minilimo-dynamic-coupled-strong-0", dh) do vtk
+VTKGridFile("minilimo-coupled-dynamic-strong-0", dh) do vtk
     write_solution(vtk, dh, u)
     Ferrite.write_node_data(vtk, resu, "ru")
     Ferrite.write_node_data(vtk, resθ, "rθ")
@@ -245,7 +245,7 @@ un = zeros(N_dof)
 #                 end
 #             end
 #             d, G3 = director_field(dh, scv, u)
-#             VTKGridFile("minilimo-dynamic-coupled-strong-$(vtk_step[])", dh) do vtk
+#             VTKGridFile("minilimo-coupled-dynamic-strong-$(vtk_step[])", dh) do vtk
 #                 write_solution(vtk, dh, u)
 #                 Ferrite.write_node_data(vtk, resu, "ru")
 #                 Ferrite.write_node_data(vtk, resθ, "rθ")
@@ -272,6 +272,7 @@ un .= load("limo_dynamic_coupled_u0.jld2")["u"]
 # Freeze the fully-morphed edge configuration (t·5 ≥ T_morph → ramp = 1) for the coupled
 # phase; the Dirichlet morph is held constant from here on (u, v, a carried forward).
 Ferrite.update!(ch, T_sim * 5)
+u .= un
 apply!(u, ch)
 
 # actuation waveform (normalized to [0,1])
@@ -289,13 +290,13 @@ wk = (; Ra, Rp, Rv, Ca, Cv, Pscale = p_max)
 
 # coupling controls
 tol_cpl  = 1e-4
-max_iter = 20
-dt_cpl   = 0.01   # doubles as the HHT-α time step in the coupled phase
+max_iter = 50
+dt_cpl   = 0.001   # doubles as the HHT-α time step in the coupled phase
 T_beat   = 4.0    # total coupled duration [s]
 
 # storages
 vols = Float64[]; pres = Float64[]; pact = Float64[]
-paos = Float64[]; pvns = Float64[]
+paos = Float64[]; pvns = Float64[]; tsav = Float64[]
 
 bufs_cpl = (; K_int, r_int, K_plv, F_plv, K_pact, F_pact, K_plvpact, F_plvpact, M, K_eff,
               res, rhs, v1, v2, dVdu, a_new, v_new, Ma, Mv, F_lu, free, g_old,
@@ -381,32 +382,50 @@ end
 V_LV0 = -2 * compute_volume(dh, scv, u; cellset=Plv_srf)   # m³
 println("Initial volume of the device: ", round(V_LV0 * m3_to_ml; digits=4), " ml")
 
-println("\nPHASE 2 — monolithic strong 3D-0D coupling (dt_cpl=$(dt_cpl) s)")
-println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters")
+println("\nPHASE 2 — monolithic strong 3D-0D coupling (adaptive Δt, Δt₀=$(dt_cpl) s)")
+println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters |    Δt [s]")
 
-# initial 0D state (Pa): filled ventricle, arterial at 80 mmHg
-@time let V_LVₙ = V_LV0, Paₙ = Pa, Pvₙ = Pv, t_cpl = 0.0, Plv=p_max, Pa = 80.0 / Pa2mmHg, Pv = p_max
+# initial 0D state [Pa]: filled ventricle (Plv=Pv=p_max), arterial at 80 mmHg
+Pa0 = 80.0 / Pa2mmHg
+Pv0 = p_max
+# Adaptive Δt driven by 3D-0D coupling convergence (same policy as the Phase-1 morph):
+# a converged monolithic step commits and grows Δt by 1.2× (capped at Δt_max); a failed
+# step is discarded (structural + 0D history untouched), Δt halved, and the step retried.
+@time let V_LVₙ = V_LV0, Paₙ = Pa0, Pvₙ = Pv0, t_cpl = 0.0, Plv = p_max, Pa = Pa0, Pv = Pv0
     step = 0
-    while t_cpl < T_beat - dt_cpl / 2
-        step += 1
-        t_cpl += dt_cpl
+    Δt_cur = dt_cpl
+    while t_cpl < T_beat - Δt_min
+        t_new  = min(t_cpl + Δt_cur, T_beat)
+        Δt_cur = t_new - t_cpl   # clip the final step to land exactly on T_beat
 
-        Pact_mmHg = 200 * ϕᵢ(t_cpl; tC=0.1, tR=0.4, TC=0.3, TR=0.3)
+        Pact_mmHg = 600 * ϕᵢ(t_new; tC=0.1, tR=0.4, TC=0.3, TR=0.3)
         Pact = Pact_mmHg / Pa2mmHg
 
-        # HHT-α predictors (Δt = dt_cpl), morph BC frozen.
-        @. ũ = u + dt_cpl * v + (dt_cpl^2 * (0.5 - β_hht)) * a
-        @. ṽ = v + (dt_cpl * (1 - γ_hht)) * a
+        # HHT-α predictors (Δt = Δt_cur), morph BC frozen.
+        @. ũ = u + Δt_cur * v + (Δt_cur^2 * (0.5 - β_hht)) * a
+        @. ṽ = v + (Δt_cur * (1 - γ_hht)) * a
         u_new .= ũ
         apply!(u_new, ch)
 
-        Plv, Pa, Pv, n_iter, converged, V₃D = solve_coupled_strong_step!(
-            u_new, ũ, ṽ, Plv, Pa, Pv, V_LVₙ, Paₙ, Pvₙ, Pact, dt_cpl, dh, scv, mat, ch,
+        # Trial pressures start from the last committed 0D state; only adopted on success.
+        Plv_t, Pa_t, Pv_t, n_iter, converged, V₃D = solve_coupled_strong_step!(
+            u_new, ũ, ṽ, Plv, Pa, Pv, V_LVₙ, Paₙ, Pvₙ, Pact, Δt_cur, dh, scv, mat, ch,
             Plv_srf, Pact_srf, PlvPact_srf, wk, bufs_cpl; max_iter=max_iter, tol=tol_cpl, verbose=false)
 
+        if !converged
+            Δt_cur /= 2
+            Δt_cur < Δt_min && error("minimum Δt reached at t=$(round(t_cpl, digits=4)) s")
+            @warn @sprintf("coupling did not converge at t=%.4f s — retrying with Δt=%.3e s", t_new, Δt_cur)
+            continue   # discard trial state (u, v, a, 0D history all unchanged), retry
+        end
+
+        step += 1
+        t_cpl = t_new
+        Plv, Pa, Pv = Plv_t, Pa_t, Pv_t
+
         # commit dynamic structural state + advance 0D history
-        @. a = (u_new - ũ) / (β_hht * dt_cpl^2)
-        @. v = ṽ + (dt_cpl * γ_hht) * a
+        @. a = (u_new - ũ) / (β_hht * Δt_cur^2)
+        @. v = ṽ + (Δt_cur * γ_hht) * a
         mul!(Mv, M, v)
         @. g_old = α_damp * Mv + r_int - (Plv * F_plv + Pact * F_pact - Pact * F_plvpact)
         u .= u_new
@@ -421,7 +440,7 @@ println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters")
             end
         end
         d, G3 = director_field(dh, scv, u)
-        VTKGridFile("minilimo-dynamic-coupled-strong-$(vtk_step[])", dh) do vtk
+        VTKGridFile("minilimo-coupled-dynamic-strong-$(vtk_step[])", dh) do vtk
             write_solution(vtk, dh, u)
             Ferrite.write_node_data(vtk, resu, "ru")
             Ferrite.write_node_data(vtk, resθ, "rθ")
@@ -430,25 +449,26 @@ println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters")
             for ID in 1:3; color(vtk, grid, "SRF_$ID"); end
             pvd[T_sim + t_cpl] = vtk
         end
-        @printf("  %9.4f | %11.4f | %14.4f | %14.4f | %d\n",
-                t_cpl, Plv * Pa2mmHg, 2V₃D * m3_to_ml, Pact_mmHg, n_iter)
-
-        !converged && (@warn "coupling step $step (t=$(t_cpl)) did not converge"; break)
+        @printf("  %9.4f | %11.4f | %14.4f | %14.4f | %5d | %.3e\n",
+                t_cpl, Plv * Pa2mmHg, 2V₃D * m3_to_ml, Pact_mmHg, n_iter, Δt_cur)
 
+        push!(tsav, t_cpl)             # actual (non-uniform) sample time [s]
         push!(vols, 2V₃D * m3_to_ml)   # full LV volume [ml]
         push!(pres, Plv * Pa2mmHg)     # LV pressure [mmHg]
         push!(pact, Pact_mmHg)
         push!(paos, Pa * Pa2mmHg)
         push!(pvns, Pv * Pa2mmHg)
+
+        Δt_cur = min(Δt_cur * 1.2, Δt_max)   # grow after a converged step
     end
 end
 close(pvd)
 
 using Plots
-times = collect(dt_cpl:dt_cpl:dt_cpl*length(pres))
+times = tsav   # adaptive Δt → non-uniform sample times
 p1 = plot(times, [vols, pres, pact, paos, pvns], xlabel="Time [s]",
           label=["Vlv" "Plv" "Pact" "Pao" "Pv"], lw=2, legend=:right)
 p2 = plot(vols, pres, label=:none, xlim=extrema(vols).+(-10,10), ylims=(0, 100),
           xlabel="Volume [ml]", ylabel="Pressure [mmHg]", lw=2, linez=times./maximum(times))
 plot(p1, p2)
-# savefig("minilimo-dynamic-coupled-strong-N$Np.png")
+# savefig("minilimo-coupled-dynamic-strong-N$Np.png")

…nilimo/limo_dynamic_coupled_transient.jl minilimo/limo_coupled_dynamic_weak.jlminilimo/limo_dynamic_coupled_transient.jl renamed to minilimo/limo_coupled_dynamic_weak.jl minilimo/limo_dynamic_coupled_transient.jl renamed to minilimo/limo_coupled_dynamic_weak.jl

@@ -4,7 +4,7 @@ import OrdinaryDiffEq as ODE
 
 # Fully-dynamic (HHT-α throughout) morphing + 3D-0D Windkessel-coupled beat of the
 # miniLIMO, on the rectangular multi-surface mesh built by `make_minilimo_grid`.  This is
-# the transient counterpart of `limo_dynamic_coupled.jl`: there the coupled phase is a
+# the transient counterpart of `limo_coupled_static_weak.jl`: there the coupled phase is a
 # quasi-static Schur solve; here inertia and damping are retained in the coupled phase too.
 #
 #   PHASE 1 — dynamic morph  (t ∈ [0, T_sim], HHT-α)
@@ -186,12 +186,12 @@ u = zeros(N_dof); apply!(u, ch)
 v = zeros(N_dof)
 a = zeros(N_dof)
 
-pvd = paraview_collection("minilimo-dynamic-coupled-transient")
+pvd = paraview_collection("minilimo-coupled-dynamic-weak")
 vtk_step = Ref(0)
 resu = zeros(3, getnnodes(dh.grid))
 resθ = zeros(2, getnnodes(dh.grid))
 d, G3 = director_field(dh, scv, u)
-VTKGridFile("minilimo-dynamic-coupled-transient-0", dh) do vtk
+VTKGridFile("minilimo-coupled-dynamic-weak-0", dh) do vtk
     write_solution(vtk, dh, u)
     Ferrite.write_node_data(vtk, resu, "ru")
     Ferrite.write_node_data(vtk, resθ, "rθ")
@@ -237,7 +237,7 @@ let t = 0.0; step = 0; Δt_cur = Δt; p = 0.0
                 end
             end
             d, G3 = director_field(dh, scv, u)
-            VTKGridFile("minilimo-dynamic-coupled-transient-$(vtk_step[])", dh) do vtk
+            VTKGridFile("minilimo-coupled-dynamic-weak-$(vtk_step[])", dh) do vtk
                 write_solution(vtk, dh, u)
                 Ferrite.write_node_data(vtk, resu, "ru")
                 Ferrite.write_node_data(vtk, resθ, "rθ")
@@ -308,7 +308,7 @@ dt_cpl   = 0.01   # doubles as the HHT-α time step in the coupled phase
 
 # storages
 vols = Float64[]; pres = Float64[]; pact = Float64[]
-paos = Float64[]; pvns = Float64[]; vtarget = Float64[]
+paos = Float64[]; pvns = Float64[]; vtarget = Float64[]; tsav = Float64[]
 
 bufs_cpl = (; K_int, r_int, K_plv, F_plv, K_pact, F_pact, K_plvpact, F_plvpact, M, K_eff,
               res, rhs, v1, v2, dVdu, a_new, v_new, Ma, Mv, F_lu, free, g_old,
@@ -366,42 +366,70 @@ end
 
 dt_cpl   = 0.005
 
-println("\nPHASE 2 — dynamic HHT-α 3D-0D coupling (dt_cpl=$(dt_cpl) s)")
-println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters")
+println("\nPHASE 2 — dynamic HHT-α 3D-0D coupling (adaptive Δt, Δt₀=$(dt_cpl) s)")
+println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters |    Δt [s]")
 
+# Adaptive Δt driven by 3D-0D coupling convergence (same policy as the Phase-1 morph):
+# a converged bordered-Newton step commits and grows Δt by 1.2× (capped at Δt_max); a failed
+# step is discarded — the 0D integrator is rolled back with `reinit!`, the structural state
+# (u, v, a) is left untouched, Δt is halved, and the step is retried.
 @time let p = p_max
     step = 0
-    while integrator.t < tspan[2] - dt_cpl / 2
-        step += 1
+    Δt_cur = dt_cpl
+    u_ode_save = zeros(length(integrator.u))
+    while integrator.t < tspan[2] - Δt_min
+        Δt_cur = min(Δt_cur, tspan[2] - integrator.t)   # clip to land exactly on tspan[2]
+
+        # snapshot the 0D state so a failed structural solve can be rolled back
+        t_ode_save = integrator.t
+        u_ode_save .= integrator.u
 
-        # advance Windkessel by dt_cpl (Plv = integrator.u[4] held fixed).
-        ODE.step!(integrator, dt_cpl, true)
+        # advance Windkessel by Δt_cur (Plv = integrator.u[4] held fixed).
+        ODE.step!(integrator, Δt_cur, true)
 
         # target half-model volume (m³) from the full-LV volume the ODE tracks.
         V_target = 0.5 * integrator.u[1] / m3_to_ml
-        push!(vtarget, integrator.u[1])
 
         # actuator pressure at this time [mmHg] → Pa
         Pact_mmHg = 200 * ϕᵢ(integrator.t; tC=0.1, tR=0.4, TC=0.3, TR=0.3)
         Pact = Pact_mmHg / Pa2mmHg
 
-        # HHT-α predictors for this step (Δt = dt_cpl), morph BC frozen.
-        @. ũ = u + dt_cpl * v + (dt_cpl^2 * (0.5 - β_hht)) * a
-        @. ṽ = v + (dt_cpl * (1 - γ_hht)) * a
+        # HHT-α predictors for this step (Δt = Δt_cur), morph BC frozen.
+        @. ũ = u + Δt_cur * v + (Δt_cur^2 * (0.5 - β_hht)) * a
+        @. ṽ = v + (Δt_cur * (1 - γ_hht)) * a
         u_new .= ũ
         apply!(u_new, ch)
 
-        p, n_iter, converged, V₃D = solve_coupled_dyn_step!(u_new, ũ, ṽ, p, Pact, V_target, dt_cpl,
+        # trial Plv starts from the last committed value; only adopted on success.
+        p_t, n_iter, converged, V₃D = solve_coupled_dyn_step!(u_new, ũ, ṽ, p, Pact, V_target, Δt_cur,
                                         dh, scv, mat, ch, Plv_srf, Pact_srf, PlvPact_srf, bufs_cpl;
                                         max_iter=max_iter, tol=tol_cpl, verbose=false)
 
+        if !converged
+            # roll the 0D integrator back to the pre-step state and retry with a smaller Δt;
+            # u, v, a and p are all left at their last committed values.
+            ODE.reinit!(integrator, u_ode_save; t0=t_ode_save)
+            @warn @sprintf("coupling did not converge for step to t=%.4f s — halving Δt %.3e → %.3e",
+                           t_ode_save + Δt_cur, Δt_cur, Δt_cur / 2)
+            Δt_cur /= 2
+            Δt_cur < Δt_min && error("minimum Δt reached at t=$(round(t_ode_save, digits=4)) s")
+            continue
+        end
+
+        step += 1
+        p = p_t
+
         # commit dynamic state (velocity/accel updates + HHT history g_old).
-        @. a = (u_new - ũ) / (β_hht * dt_cpl^2)
-        @. v = ṽ + (dt_cpl * γ_hht) * a
+        @. a = (u_new - ũ) / (β_hht * Δt_cur^2)
+        @. v = ṽ + (Δt_cur * γ_hht) * a
         mul!(Mv, M, v)
         @. g_old = α_damp * Mv + r_int - (p * F_plv + Pact * F_pact - Pact * F_plvpact)
         u .= u_new
 
+        # feed the converged LV pressure back into the ODE state.
+        integrator.u[4] = p * Pa2mmHg
+        ODE.u_modified!(integrator, true)
+
         if step%50 == 0
             vtk_step[] += 1
             for cell in CellIterator(dh)
@@ -412,7 +440,7 @@ println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters")
                 end
             end
             d, G3 = director_field(dh, scv, u)
-            VTKGridFile("minilimo-dynamic-coupled-transient-$(vtk_step[])", dh) do vtk
+            VTKGridFile("minilimo-coupled-dynamic-weak-$(vtk_step[])", dh) do vtk
                 write_solution(vtk, dh, u)
                 Ferrite.write_node_data(vtk, resu, "ru")
                 Ferrite.write_node_data(vtk, resθ, "rθ")
@@ -422,29 +450,27 @@ println("      t [s] |  p [mmHg]   |  Vlv_full [ml]  |  Pact [mmHg]  | iters")
                 pvd[T_sim + integrator.t] = vtk
             end
         end
-        @printf("  %9.4f | %11.4f | %14.4f | %14.4f | %d\n",
-                integrator.t, p * Pa2mmHg, 2V₃D * m3_to_ml, Pact_mmHg, n_iter)
-
-        !converged && (@warn "coupling step $step (t=$(integrator.t)) did not converge"; break)
-
-        # feed the converged LV pressure back into the ODE state.
-        integrator.u[4] = p * Pa2mmHg
-        ODE.u_modified!(integrator, true)
+        @printf("  %9.4f | %11.4f | %14.4f | %14.4f | %5d | %.3e\n",
+                integrator.t, p * Pa2mmHg, 2V₃D * m3_to_ml, Pact_mmHg, n_iter, Δt_cur)
 
+        push!(tsav, integrator.t)      # actual (non-uniform) sample time [s]
+        push!(vtarget, integrator.u[1])
         push!(vols, 2V₃D * m3_to_ml)   # full LV volume [ml]
         push!(pres, p * Pa2mmHg)       # LV pressure [mmHg]
         push!(pact, Pact_mmHg)
         push!(paos, integrator.u[2])
         push!(pvns, integrator.u[3])
+
+        Δt_cur = min(Δt_cur * 1.2, Δt_max)   # grow after a converged step
     end
 end
 close(pvd)
 
 using Plots
-times = collect(0:dt_cpl:integrator.t)[1:length(pres)]
+times = tsav   # adaptive Δt → non-uniform sample times
 p1 = plot(times, [vols, pres, pact, paos, pvns], xlabel="Time [s]",
           label=["Vlv" "Plv" "Pact" "Pao" "Pv"], lw=2, legend=:right)
 p2 = plot(vols, pres, label=:none, xlim=extrema(vols).+(-10,10), ylims=(0, 100),
           xlabel="Volume [ml]", ylabel="Pressure [mmHg]", lw=2, linez=times./maximum(times))
 plot(p1, p2)
-# savefig("minilimo-dynamic-coupled-transient-N$Np.png")
+# savefig("minilimo-coupled-dynamic-weak-N$Np.png")