diff --git a/Project.toml b/Project.toml index 9319887..2d82b23 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "openBF" uuid = "e815b1a4-10eb-11ea-25f1-272ff651e618" authors = ["alessandro "] -version = "2.6.0" +version = "2.6.1" [deps] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/src/anastomosis.jl b/src/anastomosis.jl index 82b2a45..84b79b8 100755 --- a/src/anastomosis.jl +++ b/src/anastomosis.jl @@ -26,7 +26,7 @@ getUan(v1::Vessel, v2::Vessel, v3::Vessel) = SVector{6, Float64}( function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel) k = (sqrt(1.5*v1.gamma[end]), sqrt(1.5*v2.gamma[end]), sqrt(1.5*v3.gamma[1])) U = getUan(v1, v2, v3) - W = (U[1] + 4k[1]* U[4], U[2] + 4k[2] * U[5], U[3] + 4k[3] * U[6]) + W = (U[1] + 4k[1]* U[4], U[2] + 4k[2] * U[5], U[3] - 4k[3] * U[6]) J = getJan(v1, v2, v3, U, k) F = getFan(v1, v2, v3, U, k, W) @@ -56,7 +56,7 @@ function getFan(v1::Vessel, v2::Vessel, v3::Vessel, U, k, W) v2.beta[end] * (U[5]^2 / sqrt(v2.A0[end]) - 1) - (v3.beta[1] * ((U[6]^2) / sqrt(v3.A0[1]) - 1))) end -function getJan(v1::Vessel, v2::Vessel, v3::Vessel, U::SArray, k::SArray) +function getJan(v1::Vessel, v2::Vessel, v3::Vessel, U, k) J::Array{Float64, 2} = zeros(Float64, 6, 6) J[1, 1] = 1.0 @@ -89,5 +89,5 @@ function NRan(U, W, J, F, k, v1, v2, v3) F = getFan(v1, v2, v3, U, k, W) J = getJan(v1, v2, v3, U, k) end - J + U end diff --git a/src/bifurcations.jl b/src/bifurcations.jl index 75763ca..ae3501d 100755 --- a/src/bifurcations.jl +++ b/src/bifurcations.jl @@ -24,7 +24,6 @@ function getJbif(v1::Vessel, v2::Vessel, v3::Vessel, U, k) J[1, 1] = 1.0 J[2, 2] = 1.0 J[3, 3] = 1.0 - J[4, 4] = 1.0 J[1, 4] = 4k[1] J[2, 5] = -4k[2] diff --git a/src/solver.jl b/src/solver.jl index 24ad28a..2773c4c 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -22,52 +22,82 @@ 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])) - maxspeed = ifelse(speed>maxspeed, speed, maxspeed) + maxspeed = max(maxspeed, speed) end Δt = v.dx / maxspeed - minΔt = ifelse(Δt < minΔt, Δt, minΔt) + minΔt = min(minΔt, Δt) end minΔt*n.Ccfl end +function skipthis(n::Network, s, t) + if n.vessels[(s, t)].solved + return true + end + + if s == 1 + return false + end + + indeg::Int64 = Graphs.indegree(n.graph, s) + + if indeg == 1 + parent_src = first(Graphs.inneighbors(n.graph, s)) + return ~n.vessels[parent_src, s].solved + end + + if indeg == 2 + parent_src = Graphs.inneighbors(n.graph, s) + return ~(n.vessels[first(parent_src), s].solved && n.vessels[last(parent_src), s].solved) + end + + # TODO: this should error + + return false +end + function solve!(n::Network, dt::Float64, current_time::Float64)::Nothing - for edge in n.edges - s = Graphs.src(edge) - t = Graphs.dst(edge) - - # inlet - s == 1 && inbc!(n.vessels[(s, t)], current_time, dt, n.heart) - - # TODO: multiple inlets - - # vessel - muscl!(n.vessels[(s, t)], dt, n.blood) - - # downstream - outdeg::Int64 = Graphs.outdegree(n.graph, t) - if outdeg == 0 # outlet - outbc!(n.vessels[(s, t)], dt, n.blood.rho) - elseif outdeg == 1 - indeg::Int64 = Graphs.indegree(n.graph, t) - d::Int64 = first(Graphs.outneighbors(n.graph, t)) - if indeg == 1 # conjunction - join_vessels!(n.vessels[(s, t)], n.vessels[(t, d)], n.blood.rho) + for (k, v) in n.vessels + v.solved = false + end + + while ~all(v.solved for v in values(n.vessels)) + for edge in n.edges + s = Graphs.src(edge) + t = Graphs.dst(edge) - # TODO: test ANASTOMOSIS!!! - elseif indeg == 2 # anastomosis - ps::Vector{Int64} = Graphs.inneighbors(n.graph, t) - if t == max(ps[1], ps[2]) + skipthis(n, s, t) && continue + + # inlet + s == 1 && inbc!(n.vessels[(s, t)], current_time, dt, n.heart) + + # TODO: multiple inlets + + # vessel + muscl!(n.vessels[(s, t)], dt, n.blood) + + # downstream + outdeg::Int64 = Graphs.outdegree(n.graph, t) + if outdeg == 0 # outlet + outbc!(n.vessels[(s, t)], dt, n.blood.rho) + elseif outdeg == 1 + indeg::Int64 = Graphs.indegree(n.graph, t) + d::Int64 = first(Graphs.outneighbors(n.graph, t)) + if indeg == 1 # conjunction + join_vessels!(n.vessels[(s, t)], n.vessels[(t, d)], n.blood.rho) + elseif indeg == 2 # anastomosis + ps::Vector{Int64} = Graphs.inneighbors(n.graph, t) solveAnastomosis( n.vessels[(ps[1], t)], n.vessels[(ps[2], t)], n.vessels[(t, d)], ) end + elseif outdeg == 2 # bifurcation + 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 - - elseif outdeg == 2 # bifurcation - ds::Vector{Int64} = Graphs.outneighbors(n.graph, t) - join_vessels!(n.vessels[(s, t)], n.vessels[t, ds[1]], n.vessels[t, ds[2]]) + n.vessels[(s, t)].solved = true end end end diff --git a/src/vessel.jl b/src/vessel.jl index ddfd077..86720b6 100644 --- a/src/vessel.jl +++ b/src/vessel.jl @@ -28,6 +28,7 @@ function Blood(config::Dict) end mutable struct Vessel + solved::Bool label::SubString{String} tosave::Bool @@ -264,7 +265,10 @@ function Vessel(config::Dict{Any,Any}, b::Blood, jump::Int64, tokeep::Vector{Str waveforms[q] = zeros(Float64, jump, 6) end + solved = false + Vessel( + solved, vessel_name, tosave, sn,