Optimise superbee! (#66)

* fastmath and remove sign * fix for a,b=0 * SArray * dt * SArray in anastomosis * 2.4.0
INSIGNEO · Apr 9, 2024 · e4535ad · e4535ad
1 parent efa23d8
commit e4535ad
Show file tree

Hide file tree

Showing 10 changed files with 97 additions and 107 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "openBF"
 uuid = "e815b1a4-10eb-11ea-25f1-272ff651e618"
 authors = ["alessandro <[email protected]>"]
-version = "2.3.0"
+version = "2.4.0"
 
 [deps]
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"

diff --git a/src/anastomosis.jl b/src/anastomosis.jl
@@ -17,7 +17,7 @@ limitations under the License.
 function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
 
     #Unknowns vector
-    U = [
+    U = @SArray [
         v1.u[end],
         v2.u[end],
         v3.u[1],
@@ -27,7 +27,7 @@ function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
     ]
 
     #Parameters vector
-    k = sqrt.(1.5 .* (v1.gamma[end], v2.gamma[end], v3.gamma[1]))
+    k = @SArray [sqrt(1.5*v1.gamma[end]), sqrt(1.5*v2.gamma[end]), sqrt(1.5*v3.gamma[1])]
     W = calculateWstarAn(U, k)
     J = calculateJacobianAn(v1, v2, v3, U, k)
     F = calculateFofUAn(v1, v2, v3, U, k, W)
@@ -40,12 +40,6 @@ function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
         dU = J \ (-F)
         U_new = U + dU
 
-        if any(isnan(dot(F, F)))
-            println(F)
-            @printf "error at anastomosis with vessels %s, %s, and %s \n" v1.label v2.label v3.label
-            break
-        end
-
         u_ok = 0
         f_ok = 0
         for i = 1:length(dU)
@@ -77,27 +71,18 @@ function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
     v1.Q[end] = v1.u[end] * v1.A[end]
     v2.Q[end] = v2.u[end] * v2.A[end]
     v3.Q[1] = v3.u[1] * v3.A[1]
-
-    v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
-    v2.P[end] = pressure(v2.A[end], v2.A0[end], v2.beta[end], v2.Pext)
-    v3.P[1] = pressure(v3.A[1], v3.A0[1], v3.beta[1], v3.Pext)
-
-    v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
-    v2.c[end] = waveSpeed(v2.A[end], v2.gamma[end])
-    v3.c[1] = waveSpeed(v3.A[1], v3.gamma[1])
-
 end
 
-function calculateWstarAn(U::Array, k::Tuple)
+function calculateWstarAn(U::SArray, k::SArray)
 
     W1 = U[1] + 4 * k[1] * U[4]
     W2 = U[2] + 4 * k[2] * U[5]
     W3 = U[3] - 4 * k[3] * U[6]
 
-    return [W1, W2, W3]
+    return @SArray [W1, W2, W3]
 end
 
-function calculateFofUAn(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tuple, W::Array)
+function calculateFofUAn(v1::Vessel, v2::Vessel, v3::Vessel, U::SArray, k::SArray, W::SArray)
 
     f1 = U[1] + 4 * k[1] * U[4] - W[1]
 
@@ -117,12 +102,12 @@ function calculateFofUAn(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tuple,
         v2.beta[end] * (U[5]^2 / sqrt(v2.A0[end]) - 1) -
         (v3.beta[1] * ((U[6]^2) / sqrt(v3.A0[1]) - 1))
 
-    return [f1, f2, f3, f4, f5, f6]
+    return @SArray [f1, f2, f3, f4, f5, f6]
 end
 
-function calculateJacobianAn(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tuple)
-
-    J = eye(6)
+function calculateJacobianAn(v1::Vessel, v2::Vessel, v3::Vessel, U::SArray, k::SArray)
+    J = @MArray zeros(6, 6) 
+    J .+= I(6)
 
     J[1, 4] = 4 * k[1]
     J[2, 5] = 4 * k[2]
@@ -142,5 +127,5 @@ function calculateJacobianAn(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tu
     J[6, 5] = 2 * v2.beta[end] * U[5] / sqrt(v2.A0[end])
     J[6, 6] = -2 * v3.beta[1] * U[6] / sqrt(v3.A0[1])
 
-    return J
+    J
 end
diff --git a/src/bifurcations.jl b/src/bifurcations.jl
@@ -16,7 +16,7 @@ limitations under the License.
 
 function join_vessels!(v1::Vessel, v2::Vessel, v3::Vessel)
     #Unknowns vector
-    U = [
+    U = @SArray [
         v1.u[end],
         v2.u[1],
         v3.u[1],
@@ -26,7 +26,7 @@ function join_vessels!(v1::Vessel, v2::Vessel, v3::Vessel)
     ]
 
     #Parameters vector
-    k = sqrt.(1.5 .* (v1.gamma[end], v2.gamma[1], v3.gamma[1]))
+    k = @SArray [sqrt(1.5*v1.gamma[end]), sqrt(1.5*v2.gamma[1]), sqrt(1.5*v3.gamma[1])]
     W = w_star_bifurcation(U, k)
     J = jacobian_bifurcation(v1, v2, v3, U, k)
     F = f_bifurcation(v1, v2, v3, U, k, W)
@@ -39,12 +39,6 @@ function join_vessels!(v1::Vessel, v2::Vessel, v3::Vessel)
         dU = J \ (-F)
         U_new = U + dU
 
-        if any(isnan(dot(F, F)))
-            println(F)
-            @printf "error at bifurcation with vessels %s, %s, and %s \n" v1.label v2.label v3.label
-            break
-        end
-
         u_ok = 0
         f_ok = 0
         for i = 1:length(dU)
@@ -76,27 +70,19 @@ function join_vessels!(v1::Vessel, v2::Vessel, v3::Vessel)
     v1.Q[end] = v1.u[end] * v1.A[end]
     v2.Q[1] = v2.u[1] * v2.A[1]
     v3.Q[1] = v3.u[1] * v3.A[1]
-
-    v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
-    v2.P[1] = pressure(v2.A[1], v2.A0[1], v2.beta[1], v2.Pext)
-    v3.P[1] = pressure(v3.A[1], v3.A0[1], v3.beta[1], v3.Pext)
-
-    v1.c[end] = wave_speed(v1.A[end], v1.gamma[end])
-    v2.c[1] = wave_speed(v2.A[1], v2.gamma[1])
-    v3.c[1] = wave_speed(v3.A[1], v3.gamma[1])
 end
 
 
-function w_star_bifurcation(U::Array, k::Tuple)
+function w_star_bifurcation(U::SArray, k::SArray)
     W1 = U[1] + 4 * k[1] * U[4]
     W2 = U[2] - 4 * k[2] * U[5]
     W3 = U[3] - 4 * k[3] * U[6]
 
-    [W1, W2, W3]
+    @SArray [W1, W2, W3]
 end
 
 
-function f_bifurcation(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tuple, W::Array)
+function f_bifurcation(v1::Vessel, v2::Vessel, v3::Vessel, U::SArray, k::SArray, W::SArray)
     f1 = U[1] + 4 * k[1] * U[4] - W[1]
     f2 = U[2] - 4 * k[2] * U[5] - W[2]
     f3 = U[3] - 4 * k[3] * U[6] - W[3]
@@ -109,12 +95,13 @@ function f_bifurcation(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tuple, W
     f6 =
         v1.beta[end] * (U[4] * U[4] / sqrt(v1.A0[end]) - 1) -
         (v3.beta[1] * (U[6] * U[6] / sqrt(v3.A0[1]) - 1))
-    [f1, f2, f3, f4, f5, f6]
+    @SArray [f1, f2, f3, f4, f5, f6]
 end
 
 
-function jacobian_bifurcation(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Tuple)
-    J = zeros(6, 6) + I(6)
+function jacobian_bifurcation(v1::Vessel, v2::Vessel, v3::Vessel, U::SArray, k::SArray)
+    J = @MArray zeros(6, 6) 
+    J .+= I(6)
 
     J[1, 4] = 4 * k[1]
     J[2, 5] = -4 * k[2]

diff --git a/src/boundary_conditions.jl b/src/boundary_conditions.jl
@@ -35,26 +35,28 @@ function inlet_from_data(t::Float64, h::Heart)
 end
 
 function riemann_invariants(i::Int64, v::Vessel)
-    W1 = v.u[i] - 4 * v.c[i]
-    W2 = v.u[i] + 4 * v.c[i]
-    W1, W2
+    c = wave_speed(v.A[i], v.gamma[i])
+    # W1 = v.u[i] - 4c
+    # W2 = v.u[i] + 4c
+    # W1, W2
+    @SArray [v.u[i] - 4c, v.u[i] + 4c]
 end
 
 function inv_riemann_invariants(W1::Float64, W2::Float64)
-    u = 0.5 * (W1 + W2)
-    c = (W2 - W1) * 0.125
-    u, c
+    0.5 * (W1 + W2) # u
+    # c = (W2 - W1) * 0.125
+    # u, c
 end
 
 function inlet_compatibility(dt::Float64, v::Vessel)
 
     W11, W21 = riemann_invariants(1, v)
     W12, W22 = riemann_invariants(2, v)
 
-    W11 += (W12 - W11) * (v.c[1] - v.u[1]) * dt / v.dx
+    W11 += (W12 - W11) * (wave_speed(v.A[1], v.gamma[1]) - v.u[1]) * dt / v.dx
     W21 = 2 * v.Q[1] / v.A[1] - W11
 
-    v.u[1], v.c[1] = inv_riemann_invariants(W11, W21)
+    v.u[1] = inv_riemann_invariants(W11, W21)
 
     v.A[1] = v.Q[1] / v.u[1]
     v.P[1] = pressure(v.A[1], v.A0[1], v.beta[1], v.Pext)
@@ -77,10 +79,10 @@ function outlet_compatibility(dt::Float64, v::Vessel)
     W1M1, W2M1 = riemann_invariants(v.M - 1, v)
     W1M, W2M = riemann_invariants(v.M, v)
 
-    W2M += (W2M1 - W2M) * (v.u[end] + v.c[end]) * dt / v.dx
+    W2M += (W2M1 - W2M) * (v.u[end] + wave_speed(v.A[end], v.gamma[end])) * dt / v.dx
     W1M = v.W1M0 - v.Rt * (W2M - v.W2M0)
 
-    v.u[end], v.c[end] = inv_riemann_invariants(W1M, W2M)
+    v.u[end] = inv_riemann_invariants(W1M, W2M)
     v.Q[end] = v.A[end] * v.u[end]
 end
 

diff --git a/src/conjunctions.jl b/src/conjunctions.jl
@@ -16,8 +16,8 @@ limitations under the License.
 
 function join_vessels!(b::Blood, v1::Vessel, v2::Vessel)
 
-    U = [v1.u[end], v2.u[1], sqrt(sqrt(v1.A[end])), sqrt(sqrt(v2.A[1]))]
-    k = sqrt.(1.5 .* (v1.gamma[end], v2.gamma[1]))
+    U = @SArray [v1.u[end], v2.u[1],  sqrt(sqrt(v1.A[end])), sqrt(sqrt(v2.A[1]))]
+    k = @SArray [sqrt(1.5*v1.gamma[end]), sqrt(1.5*v2.gamma[1])]
     W = w_star_conj(U, k)
     J = jacobian_conj(b, v1, v2, U, k)
     F = f_conj(b, v1, v2, U, k, W)
@@ -29,11 +29,6 @@ function join_vessels!(b::Blood, v1::Vessel, v2::Vessel)
         dU = J \ (-F)
         U_new = U + dU
 
-        if any(isnan(dot(F, F)))
-            println(F)
-            break
-        end
-
         u_ok = 0
         f_ok = 0
         for i = 1:length(dU)
@@ -61,25 +56,19 @@ function join_vessels!(b::Blood, v1::Vessel, v2::Vessel)
 
     v2.A[1] = U[4] * U[4] * U[4] * U[4]
     v2.Q[1] = v2.u[1] * v2.A[1]
-
-    v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
-    v2.P[1] = pressure(v2.A[1], v2.A0[1], v2.beta[1], v2.Pext)
-
-    v1.c[end] = wave_speed(v1.A[end], v1.gamma[end])
-    v2.c[1] = wave_speed(v2.A[1], v2.gamma[1])
 end
 
 
-function w_star_conj(U::Array, k::Tuple)
+function w_star_conj(U::SArray, k::SArray)
 
     W1 = U[1] + 4 * k[1] * U[3]
     W2 = U[2] - 4 * k[2] * U[4]
 
-    [W1, W2]
+    @SArray [W1, W2]
 end
 
 
-function f_conj(b::Blood, v1::Vessel, v2::Vessel, U::Array, k::Tuple, W::Array)
+function f_conj(b::Blood, v1::Vessel, v2::Vessel, U::SArray, k::SArray, W::SArray)
 
     f1 = U[1] + 4 * k[1] * U[3] - W[1]
 
@@ -91,13 +80,14 @@ function f_conj(b::Blood, v1::Vessel, v2::Vessel, U::Array, k::Tuple, W::Array)
         0.5 * b.rho * U[1] * U[1] + v1.beta[end] * (U[3] * U[3] / sqrt(v1.A0[end]) - 1) -
         (0.5 * b.rho * U[2] * U[2] + v2.beta[1] * (U[4] * U[4] / sqrt(v2.A0[1]) - 1))
 
-    [f1, f2, f3, f4]
+    @SArray [f1, f2, f3, f4]
 end
 
 
-function jacobian_conj(b::Blood, v1::Vessel, v2::Vessel, U::Array, k::Tuple)
+function jacobian_conj(b::Blood, v1::Vessel, v2::Vessel, U::SArray, k::SArray)
 
-    J = zeros(4, 4) + I(4)
+    J = @MArray zeros(4, 4)
+    J .+= I(4)
 
     J[1, 3] = 4 * k[1]
     J[2, 4] = -4 * k[2]

diff --git a/src/network.jl b/src/network.jl
@@ -31,6 +31,7 @@ struct Network
     vessels::Dict{Tuple{Int,Int},Vessel}
     blood::Blood
     heart::Heart
+    Ccfl::Float64
 end
 number_of_nodes(config::Vector{Dict{Any,Any}}) = maximum(c["tn"] for c in config)
 function Network(
@@ -48,13 +49,13 @@ function Network(
 
     vessels = Dict()
     for vessel_config in config
-        vessel = Vessel(vessel_config, blood, Ccfl, jump, tokeep)
+        vessel = Vessel(vessel_config, blood, jump, tokeep)
         add_edge!(graph, vessel.sn, vessel.tn)
         vessels[(vessel.sn, vessel.tn)] = vessel
         verbose && next!(prog)
     end
     check(graph)
-    Network(graph, vessels, blood, heart)
+    Network(graph, vessels, blood, heart, Ccfl)
 end
 
 function check(g::SimpleDiGraph)

diff --git a/src/openBF.jl b/src/openBF.jl
@@ -20,7 +20,7 @@ module openBF
 using Printf
 using LinearAlgebra
 using DelimitedFiles
-
+using StaticArrays
 using YAML
 using Glob
 using Graphs

diff --git a/src/simulation.jl b/src/simulation.jl
@@ -122,7 +122,7 @@ function run_simulation(
     stats = @timed while true
 
         # step
-        dt = calculate_Δt(network)
+        dt = calculateΔt(network)
         solve!(network, dt, current_time)
         update_ghost_cells!(network)
         verbose && next!(prog)