match inlet impedance

src/boundary_conditions.jl

@@ -62,16 +62,14 @@ function inlet_compatibility(dt::Float64, v::Vessel)
 end
 
 
-function set_outlet_bc(dt::Float64, v::Vessel)
+function set_outlet_bc(dt::Float64, v::Vessel, ρ::Float64)
     # A proximal resistance `R1` set to zero means that a reflection coefficient
     if v.R1 == 0.0
         v.P[end] = 2 * v.P[end-1] - v.P[end-2]
         outlet_compatibility(dt, v)
     else
-        wk3(dt, v)
-
+        wk3(dt, v, ρ)
     end
-
 end
 
 function outlet_compatibility(dt::Float64, v::Vessel)
@@ -87,10 +85,15 @@ function outlet_compatibility(dt::Float64, v::Vessel)
 end
 
 
-function wk3(dt::Float64, v::Vessel)
+function wk3(dt::Float64, v::Vessel, ρ::Float64)
     Pout = 0.0
     Al = v.A[end]
     ul = v.u[end]
+
+    # inlet impedance matching
+    v.R1 = ρ * wave_speed(v.A[end], v.gamma[end]) / v.A[end]
+    v.R2 = v.total_peripheral_resistance - v.R1
+
     v.Pc += dt / v.Cc * (Al * ul - (v.Pc - Pout) / v.R2)
     As = Al
 
@@ -111,13 +114,20 @@ function wk3(dt::Float64, v::Vessel)
     v.u[end] = us
 end
 
-function newtone(f::Function, df::Function, x0)
-    xn = x0 - f(x0) / df(x0)
-    if abs(xn - x0) <= 1e-5
-        return xn
-    else
-        newtone(f, df, xn)
+# function newtone(f::Function, df::Function, x0)
+#     xn = x0 - f(x0) / df(x0)
+#     if abs(xn - x0) <= 1e-5
+#         return xn
+#     else
+#         newtone(f, df, xn)
+#     end
+# end
+
+function newtone(f::Function, df::Function, xn)
+    for _=1:10
+        xn -= f(xn) / df(xn)
     end
+    xn
 end
 
 update_ghost_cells!(n::Network) = update_ghost_cells!.(values(n.vessels))

src/solver.jl

@@ -31,7 +31,7 @@ function solve!(n::Network, dt::Float64, current_time::Float64)
         muscl!(n.vessels[(s, t)], dt, n.blood)
 
         if outdegree(n.graph, t) == 0 # outlet
-            set_outlet_bc(dt, n.vessels[(s, t)])
+            set_outlet_bc(dt, n.vessels[(s, t)], n.blood.rho)
         elseif outdegree(n.graph, t) == 1
             if indegree(n.graph, t) == 1 # conjunction
                 d = first(outneighbors(n.graph, t))

src/vessel.jl

@@ -84,6 +84,7 @@ mutable struct Vessel
 
     R1::Float64
     R2::Float64
+    total_peripheral_resistance::Float64
     Cc::Float64
     Pc::Float64
 
@@ -228,6 +229,7 @@ function Vessel(config::Dict{Any,Any}, b::Blood, Ccfl::Float64)
     if R2 == 0.0
         R2 = R1 - b.rho * wave_speed(A0[end], gamma[end]) / A0[end]
     end
+    total_peripheral_resistance = R1 + R2
 
     # `Pc` is the pressure through the peripheral compliance of the three
     # elements windkessel. It is set to zero to simulate the pressure at the
@@ -299,6 +301,7 @@ function Vessel(config::Dict{Any,Any}, b::Blood, Ccfl::Float64)
         Rt,
         R1,
         R2,
+        total_peripheral_resistance,
         Cc,
         Pc,
         slope,