Skip to content

Commit

Permalink
wip
Browse files Browse the repository at this point in the history
  • Loading branch information
@alemelis
alemelis committed Nov 23, 2023
1 parent dd65c1c commit d31ce84
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 61 deletions.
102 changes: 45 additions & 57 deletions src/anastomosis.jl
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@

function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
function join_anastomosis!(v1::Vessel, v2::Vessel, v3::Vessel)

#Unknowns vector
U = [
U = SA[
v1.u[end],
v2.u[end],
v3.u[1],
Expand All @@ -12,45 +12,33 @@ function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
]

#Parameters vector
k1 = sqrt(0.5 * 3 * v1.gamma)
k2 = sqrt(0.5 * 3 * v2.gamma)
k3 = sqrt(0.5 * 3 * v3.gamma)
k = [k1, k2, k3]
k = SA[sqrt(1.5 * v1.gamma), sqrt(1.5 * v2.gamma), sqrt(1.5 * v3.gamma)]
J = jacobian_anastomosis(v1, v2, v3, U, k)

W = calculateWstarAn(U, k)
J = calculateJacobianAn(v1, v2, v3, U, k)
F = calculateFofUAn(v1, v2, v3, U, k, W)
W = w_star_anastomosis(U, k)
F = f_anastomosis(v1, v2, v3, U, k, W)

#Newton-Raphson
nr_toll_U = 1.e-5
nr_toll_F = 1.e-5

while true
dU = J \ (-F)
# U_new = U + 0.01*dU
U_new = U + dU
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)
if abs(dU[i]) <= nr_toll_U || abs(F[i]) <= nr_toll_F
u_ok += 1
f_ok += 1
end
end
u_ok = sum(abs.(dU) .<= nr_toll_U)
f_ok = sum(abs.(F) .<= nr_toll_F)

if u_ok == length(dU) || f_ok == length(dU)
U = U_new
if u_ok == 6 || f_ok == 6
break
else
U = U_new
W = calculateWstarBif(U, k)
W = w_star_anastomosis(U, k)
F = calculateFofUBif(v1, v2, v3, U, k, W)
end
end
Expand All @@ -72,65 +60,65 @@ function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
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])
v1.c[end] = wave_speed(v1.A[end], v1.gamma[end])
v2.c[end] = wave_speed(v2.A[end], v2.gamma[end])
v3.c[1] = wave_speed(v3.A[1], v3.gamma[1])

end

function calculateWstarAn(U::Array, k::Array)

W1 = U[1] + 4 * k[1] * U[4]
W2 = U[2] + 4 * k[2] * U[5]
W3 = U[3] - 4 * k[3] * U[6]
function w_star_anastomosis(U::SVector, k::SVector)
W1 = U[1] + 4k[1] * U[4]
W2 = U[2] + 4k[2] * U[5]
W3 = U[3] - 4k[3] * U[6]

return [W1, W2, W3]
@SArray [W1, W2, W3]
end

function calculateFofUAn(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Array, W::Array)

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]

function f_anastomosis(
v1::Vessel,
v2::Vessel,
v3::Vessel,
U::SVector,
k::SVector,
W::SVector,
)

f1 = U[1] + 4k[1] * U[4] - W[1]
f2 = U[2] + 4k[2] * U[5] - W[2]
f3 = U[3] - 4k[3] * U[6] - W[3]
f4 =
U[1] * (U[4] * U[4] * U[4] * U[4]) + U[2] * (U[5] * U[5] * U[5] * U[5]) -
U[3] * (U[6] * U[6] * U[6] * U[6])

f5 =
v1.beta[end] * (U[4]^2 / sqrt(v1.A0[end]) - 1) -
(v3.beta[1] * ((U[6])^2 / sqrt(v3.A0[1]) - 1))

f6 =
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]
@SArray [f1, f2, f3, f4, f5, f6]
end

function calculateJacobianAn(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Array)

J = eye(6)
function jacobian_anastomosis(v1::Vessel, v2::Vessel, v3::Vessel, U::Array, k::Array)
J = zeros(6, 6) + I(6)

J[1, 4] = 4 * k[1]
J[2, 5] = 4 * k[2]
J[3, 6] = -4 * k[3]
J[1, 4] = 4k[1]
J[2, 5] = 4k[2]
J[3, 6] = -4k[3]

J[4, 1] = U[4] * U[4] * U[4] * U[4]
J[4, 2] = U[5] * U[5] * U[5] * U[5]
J[4, 3] = -U[6] * U[6] * U[6] * U[6]
J[4, 4] = 4 * U[1] * (U[4] * U[4] * U[4])
J[4, 5] = 4 * U[2] * (U[5] * U[5] * U[5])
J[4, 6] = -4 * U[3] * (U[6] * U[6] * U[6])
J[4, 4] = 4U[1] * (U[4] * U[4] * U[4])
J[4, 5] = 4U[2] * (U[5] * U[5] * U[5])
J[4, 6] = -4U[3] * (U[6] * U[6] * U[6])

J[5, 4] = 2 * v1.beta[end] * U[4] / sqrt(v1.A0[end])
J[5, 4] = 2v1.beta[end] * U[4] / sqrt(v1.A0[end])
J[5, 5] = 0.0
J[5, 6] = -2 * v3.beta[1] * U[6] / sqrt(v3.A0[1])
J[5, 6] = -2v3.beta[1] * U[6] / sqrt(v3.A0[1])

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])
J[6, 5] = 2v2.beta[end] * U[5] / sqrt(v2.A0[end])
J[6, 6] = -2v3.beta[1] * U[6] / sqrt(v3.A0[1])

return J
SMatrix{6,6}(J)
end
2 changes: 1 addition & 1 deletion src/bifurcations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ function newton_raphson(
u_ok = sum(abs.(dU) .<= nr_toll_U)
f_ok = sum(abs.(F) .<= nr_toll_F)

if u_ok == length(dU) || f_ok == length(dU)
if u_ok == 6 || f_ok == 6
return U
else
W = w_star_bifurcation(U, k)
Expand Down
3 changes: 1 addition & 2 deletions src/conjunctions.jl
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@
# TODO: generic implementation for all junctions
function join_vessels!(b::Blood, v1::Vessel, v2::Vessel)

U = SA[v1.u[end], v2.u[1], sqrt(sqrt(v1.A[end])), sqrt(sqrt(v2.A[1]))]
Expand All @@ -25,7 +24,7 @@ function join_vessels!(b::Blood, v1::Vessel, v2::Vessel)
u_ok = sum(abs.(dU) .<= nr_toll_U)
f_ok = sum(abs.(F) .<= nr_toll_F)

if u_ok == length(dU) || f_ok == length(dU)
if u_ok == 4 || f_ok == 4
break
else
W = w_star_conj(U, k)
Expand Down
2 changes: 1 addition & 1 deletion src/solver.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ function solve!(n::Network, dt::Float64, current_time::Float64)
p1, p2 = inneighbors(n.graph, t)
if t == max(p1, p2)
d = outneighbors(grado, t)
solveAnastomosis(
join_anastomosis!(
n.vessels[(p1, t)],
n.vessels[(p2, t)],
n.vessels[(t, d)],
Expand Down

0 comments on commit d31ce84

Please sign in to comment.