Apply boundary conditions first and then solve the vessel (#71)

* apply boundary conditions first and then solve the vessel

* skip nans

Project.toml

@@ -1,7 +1,7 @@
 name = "openBF"
 uuid = "e815b1a4-10eb-11ea-25f1-272ff651e618"
 authors = ["alessandro <[email protected]>"]
-version = "2.6.1"
+version = "2.6.2"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/solver.jl

@@ -22,9 +22,11 @@ function calculateΔt(n::Network)
         maxspeed = 0.0
         for i=eachindex(v.u)
             @inbounds speed = abs(v.u[i] + wave_speed(v.A[i], v.gamma[i+1]))
+            isnan(speed) && continue
             maxspeed = max(maxspeed, speed)
         end
         Δt = v.dx / maxspeed
+        isnan(Δt) && continue
         minΔt = min(minΔt, Δt)
     end
     minΔt*n.Ccfl
@@ -73,9 +75,6 @@ function solve!(n::Network, dt::Float64, current_time::Float64)::Nothing
 
             # TODO: multiple inlets
 
-            # vessel
-            muscl!(n.vessels[(s, t)], dt, n.blood)
-
             # downstream
             outdeg::Int64 = Graphs.outdegree(n.graph, t)
             if outdeg == 0 # outlet
@@ -97,6 +96,8 @@ function solve!(n::Network, dt::Float64, current_time::Float64)::Nothing
                 ds::Vector{Int64} = Graphs.outneighbors(n.graph, t)
                 join_vessels!(n.vessels[(s, t)], n.vessels[t, ds[1]], n.vessels[t, ds[2]])
             end
+
+            muscl!(n.vessels[(s, t)], dt, n.blood)
             n.vessels[(s, t)].solved = true
         end
     end