Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Catheter aq #20

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Catheter aq #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
98bbb6d
first commit catheter branch
ibenemerito88 Dec 20, 2019
065bfbd
Ac scalar
ibenemerito88 Dec 26, 2019
6894295
add bifurcations, conjunctions
ibenemerito88 Jan 4, 2020
7d281be
modified c viscous term
ibenemerito88 Jan 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion .gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,6 @@ test/**/*_results/*
**/.DS_Store
cvs_results/*
docs/build/*
Manifest.toml
Manifest.toml
*.sh
*.swp
Copy link
Collaborator

@alemelis alemelis Jan 4, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

qui conviene aggiungere anche tutti i .last e .out così da non averli anche nella PR, insieme a tutta la cartella backupres

1 change: 1 addition & 0 deletions Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,3 +12,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
56 changes: 34 additions & 22 deletions src/anastomosis.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,19 +46,20 @@ end
Return the Jacobian for anastomosis equations.
"""
function calculateJacobianAnastomosis(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel, U, k)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

qui, ma anche nelle funzioni successive, per qualche motivo che non ricordo non ho specificato il tipo di U e k, ma andrebbe fatto. A occhio mi sembrano due Array{Float,1}

@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[4]^4) # MODIFIED THIS LINE
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

qui puoi usare U43 e calcolare g1 come

g1 = sqrt(1.0 - v1.Ac / (U43*U[4]))

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ovviamente g va definita dopo U

Copy link
Collaborator

@alemelis alemelis Jan 5, 2020 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

intendo dire che le linee da 49 a 54 vanno riorganizzate come

U43 = U[4]*U[4]*U[4]
g1 = sqrt(1.0 - v1.Ac / (U43*U[4]))
U53 = U[5]*U[5]*U[5]
g2 = sqrt(1.0 - v2.Ac/ (U53*U[5]))
U63 = U[6]*U[6]*U[6]
g3 = sqrt(1.0 - v3.Ac/U63*U[6])

@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[5]^4) # MODIFIED THIS LINE
@fastmath @inbounds g3 = sqrt(1-v3.Ac/U[6]^4) # MODIFIED THIS LINE
@fastmath @inbounds U43 = U[4]*U[4]*U[4]
@fastmath @inbounds U53 = U[5]*U[5]*U[5]
@fastmath @inbounds U63 = U[6]*U[6]*U[6]

@fastmath @inbounds J14 = 4.0*k[1]
@fastmath @inbounds J14 = 4.0*k[1]*(v1.Ac/(2*U[4]*g1) + g1) # MODIFIED THIS LINE
@fastmath @inbounds J25 = 4.0*k[2]*(v2.Ac/(2*U[5]*g2) + g2) # MODIFIED THIS LINE
@fastmath @inbounds J36 = -4.0*k[3]*(v3.Ac/(2*U[6]*g3) + g3) # MODIFIED THIS LINE

@fastmath @inbounds J25 = 4.0*k[2]

@fastmath @inbounds J36 = -4.0*k[3]

@fastmath @inbounds J41 = U[4]*U43
@fastmath @inbounds J42 = U[5]*U53
@fastmath @inbounds J43 = -U[6]*U63
@fastmath @inbounds J41 = U[4]*U43 - v1.Ac # MODIFIED THIS LINE
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

e dato che qui si ricalcola U[4]*U43, conviene definire una variabile U44 per questo valore in modo da calcolarla una volta sola

@fastmath @inbounds J42 = U[5]*U53 - v2.Ac # MODIFIED THIS LINE
@fastmath @inbounds J43 = -U[6]*U63 + v3.Ac # MODIFIED THIS LINE
@fastmath @inbounds J44 = 4.0*U[1]*U43
@fastmath @inbounds J45 = 4.0*U[2]*U53
@fastmath @inbounds J46 = -4.0*U[3]*U63
Expand All @@ -83,10 +84,17 @@ end

Return the Riemann invariants at the anastomosis node.
"""
function calculateWstarAnastomosis(U, k)
@fastmath @inbounds W1 = U[1] + 4*k[1]*U[4]
@fastmath @inbounds W2 = U[2] + 4*k[2]*U[5]
@fastmath @inbounds W3 = U[3] - 4*k[3]*U[6]
function calculateWstarAnastomosis(U, k, vessels :: Array{Vessel,1})
# MODIFIED THIS FUNCTION
v1 = vessels[1] # MODIFIED THIS LINE
v2 = vessels[2] # MODIFIED THIS LINE
v3 = vessels[3] # MODIFIED THIS LINE
@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[4]^4) # MODIFIED THIS LINE
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

in generale * é più veloce di ^, quindi puoi riscrivere come v1.Ac/(U[4]*U[4]*U[4]*U[4])

in generale se usi direttamente 1.0 invece di 1, eviti la conversione durante l'esecuzione

@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[5]^4) # MODIFIED THIS LINE
@fastmath @inbounds g3 = sqrt(1-v3.Ac/U[6]^4) # MODIFIED THIS LINE
@fastmath @inbounds W1 = U[1] + 4*k[1]*U[4]*g1 # MODIFIED THIS LINE
@fastmath @inbounds W2 = U[2] + 4*k[2]*U[5]*g2 # MODIFIED THIS LINE
@fastmath @inbounds W3 = U[3] - 4*k[3]*U[6]*g3 # MODIFIED THIS LINE

return @SArray [W1, W2, W3]
end
Expand All @@ -100,17 +108,21 @@ function calculateFAnastomosis(vessels :: Array{Vessel,1}, U, k, W)
v2 = vessels[2]
v3 = vessels[3]

@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[4]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[5]^4) # MODIFIED THIS LINE
@fastmath @inbounds g3 = sqrt(1-v3.Ac/U[6]^4) # MODIFIED THIS LINE

@fastmath @inbounds U42 = U[4]*U[4]
@fastmath @inbounds U52 = U[5]*U[5]
@fastmath @inbounds U62 = U[6]*U[6]

@fastmath @inbounds f1 = U[1] + 4*k[1]*U[4] - W[1]
@fastmath @inbounds f1 = U[1] + 4*k[1]*U[4]*g1 - W[1] # MODIFIED THIS LINE

@fastmath @inbounds f2 = U[2] + 4*k[2]*U[5] - W[2]
@fastmath @inbounds f2 = U[2] + 4*k[2]*U[5]*g2 - W[2] # MODIFIED THIS LINE

@fastmath @inbounds f3 = U[3] - 4*k[3]*U[6] - W[3]
@fastmath @inbounds f3 = U[3] - 4*k[3]*U[6]*g3 - W[3] # MODIFIED THIS LINE

@fastmath @inbounds f4 = U[1]*U42*U42 + U[2]*U52*U52 - U[3]*U62*U62
@fastmath @inbounds f4 = U[1]*(U42*U42 - v1.Ac) + U[2]*(U52*U52 - v2.Ac) - U[3]*(U62*U62 - v3.Ac) # MODIFIED THIS LINE

f5 = v1.beta[end]*(U42*v1.s_inv_A0[end] - 1.0) -
(v3.beta[1]*(U62*v3.s_inv_A0[1] - 1.0))
Expand All @@ -136,15 +148,15 @@ function updateAnastomosis(U, v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)
@fastmath @inbounds v2.A[end] = U[5]*U[5]*U[5]*U[5]
@fastmath @inbounds v3.A[1] = U[6]*U[6]*U[6]*U[6]

@fastmath @inbounds v1.Q[end] = v1.u[end]*v1.A[end]
@fastmath @inbounds v2.Q[end] = v2.u[end]*v2.A[end]
@fastmath @inbounds v3.Q[1] = v3.u[1]*v3.A[1]
@fastmath @inbounds v1.Q[end] = v1.u[end]*(v1.A[end] - v1.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v2.Q[end] = v2.u[end]*(v2.A[end] - v2.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v3.Q[1] = v3.u[1]*(v3.A[1] - v3.Ac) # MODIFIED THIS LINE

@inbounds v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
@inbounds v2.P[end] = pressure(v2.A[end], v2.A0[end], v2.beta[end], v2.Pext)
@inbounds v3.P[1] = pressure(v3.A[1], v3.A0[1], v3.beta[1], v3.Pext)

@fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
@fastmath @inbounds v2.c[end] = waveSpeed(v2.A[end], v2.gamma[end])
@fastmath @inbounds v3.c[1] = waveSpeed(v3.A[1], v3.gamma[1])
@fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end], v1.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v2.c[end] = waveSpeed(v2.A[end], v2.gamma[end], v2.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v3.c[1] = waveSpeed(v3.A[1], v3.gamma[1], v3.Ac) # MODIFIED THIS LINE
end
54 changes: 34 additions & 20 deletions src/bifurcations.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,17 +46,20 @@ end
Return the Jacobian for bifurcation equations.
"""
function calculateJacobianBifurcation(v1 :: Vessel, v2 :: Vessel, v3 :: Vessel, U, k)
@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[4]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[5]^4) # MODIFIED THIS LINE
@fastmath @inbounds g3 = sqrt(1-v3.Ac/U[6]^4) # MODIFIED THIS LINE
@fastmath @inbounds U43 = U[4]*U[4]*U[4]
@fastmath @inbounds U53 = U[5]*U[5]*U[5]
@fastmath @inbounds U63 = U[6]*U[6]*U[6]

@fastmath @inbounds J14 = 4.0*k[1]
@fastmath @inbounds J25 = -4.0*k[2]
@fastmath @inbounds J36 = -4.0*k[3]
@fastmath @inbounds J14 = 4.0*k[1]*(v1.Ac/(2*U[4]*g1) + g1) # MODIFIED THIS LINE
@fastmath @inbounds J25 = -4.0*k[2]*(v2.Ac/(2*U[5]*g2) + g2) # MODIFIED THIS LINE
@fastmath @inbounds J36 = -4.0*k[3]*(v3.Ac/(2*U[6]*g3) + g3) # MODIFIED THIS LINE

@fastmath @inbounds J41 = U[4]*U43
@fastmath @inbounds J42 = -U[5]*U53
@fastmath @inbounds J43 = -U[6]*U63
@fastmath @inbounds J41 = U[4]*U43 - v1.Ac # MODIFIED THIS LINE
@fastmath @inbounds J42 = -U[5]*U53 + v2.Ac # MODIFIED THIS LINE
@fastmath @inbounds J43 = -U[6]*U63 + v3.Ac # MODIFIED THIS LINE
@fastmath @inbounds J44 = 4.0*U[1]*U43
@fastmath @inbounds J45 = -4.0*U[2]*U53
@fastmath @inbounds J46 = -4.0*U[3]*U63
Expand All @@ -81,10 +84,17 @@ end

Return the Riemann invariants at the bifurcation node.
"""
function calculateWstarBifurcation(U, k)
@fastmath @inbounds W1 = U[1] + 4.0*k[1]*U[4]
@fastmath @inbounds W2 = U[2] - 4.0*k[2]*U[5]
@fastmath @inbounds W3 = U[3] - 4.0*k[3]*U[6]
function calculateWstarBifurcation(U, k, vessels :: Array{Vessel,1})
# MODIFIED THIS FUNCTION
v1 = vessels[1] # MODIFIED THIS LINE
v2 = vessels[2] # MODIFIED THIS LINE
v3 = vessels[3] # MODIFIED THIS LINE
@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[4]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[5]^4) # MODIFIED THIS LINE
@fastmath @inbounds g3 = sqrt(1-v3.Ac/U[6]^4) # MODIFIED THIS LINE
@fastmath @inbounds W1 = U[1] + 4.0*k[1]*U[4]*g1 # MODIFIED THIS LINE
@fastmath @inbounds W2 = U[2] - 4.0*k[2]*U[5]*g2 # MODIFIED THIS LINE
@fastmath @inbounds W3 = U[3] - 4.0*k[3]*U[6]*g3 # MODIFIED THIS LINE

return @SArray [W1, W2, W3]
end
Expand All @@ -98,17 +108,21 @@ function calculateFBifurcation(vessels :: Array{Vessel,1}, U, k, W)
v2 = vessels[2]
v3 = vessels[3]

@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[4]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[5]^4) # MODIFIED THIS LINE
@fastmath @inbounds g3 = sqrt(1-v3.Ac/U[6]^4) # MODIFIED THIS LINE

@fastmath @inbounds U42 = U[4]*U[4]
@fastmath @inbounds U52 = U[5]*U[5]
@fastmath @inbounds U62 = U[6]*U[6]

@fastmath @inbounds f1 = U[1] + 4*k[1]*U[4] - W[1]
@fastmath @inbounds f1 = U[1] + 4*k[1]*U[4]*g1 - W[1] # MODIFIED THIS LINE

@fastmath @inbounds f2 = U[2] - 4*k[2]*U[5] - W[2]
@fastmath @inbounds f2 = U[2] - 4*k[2]*U[5]*g2 - W[2] # MODIFIED THIS LINE

@fastmath @inbounds f3 = U[3] - 4*k[3]*U[6] - W[3]
@fastmath @inbounds f3 = U[3] - 4*k[3]*U[6]*g3 - W[3] # MODIFIED THIS LINE

@fastmath @inbounds f4 = U[1]*(U42*U42) - U[2]*(U52*U52) - U[3]*(U62*U62)
@fastmath @inbounds f4 = U[1]*(U42*U42 - v1.Ac) - U[2]*(U52*U52 - v2.Ac) - U[3]*(U62*U62 - v3.Ac) # MODIFIED THIS LINE

f5 = v1.beta[end]*(U42*v1.s_inv_A0[end] - 1.0) -
(v2.beta[1]*(U52*v2.s_inv_A0[1] - 1.0))
Expand All @@ -134,15 +148,15 @@ function updateBifurcation(U, v1 :: Vessel, v2 :: Vessel, v3 :: Vessel)
@fastmath @inbounds v2.A[1] = U[5]*U[5]*U[5]*U[5]
@fastmath @inbounds v3.A[1] = U[6]*U[6]*U[6]*U[6]

@fastmath @inbounds v1.Q[end] = v1.u[end]*v1.A[end]
@fastmath @inbounds v2.Q[1] = v2.u[1]*v2.A[1]
@fastmath @inbounds v3.Q[1] = v3.u[1]*v3.A[1]
@fastmath @inbounds v1.Q[end] = v1.u[end]*(v1.A[end] - v1.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v2.Q[1] = v2.u[1]*(v2.A[1] - v2.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v3.Q[1] = v3.u[1]*(v3.A[1] - v3.Ac) # MODIFIED THIS LINE

@inbounds v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
@inbounds v2.P[1] = pressure(v2.A[1], v2.A0[1], v2.beta[1], v2.Pext)
@inbounds v3.P[1] = pressure(v3.A[1], v3.A0[1], v3.beta[1], v3.Pext)

@fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
@fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1])
@fastmath @inbounds v3.c[1] = waveSpeed(v3.A[1], v3.gamma[1])
@fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end], v1.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1], v2.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v3.c[1] = waveSpeed(v3.A[1], v3.gamma[1], v3.Ac) # MODIFIED THIS LINE
end
30 changes: 15 additions & 15 deletions src/boundary_conditions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,17 +71,17 @@ function inletCompatibility(dt :: Float64, v :: Vessel, h :: Heart)
W11, W21 = riemannInvariants(1, v)
W12, W22 = riemannInvariants(2, v)

@fastmath @inbounds W11 += (W12 - W11)*(v.c[1] - v.u[1])*dt*v.invDx
@fastmath @inbounds W21 = 2.0*v.Q[1]/v.A[1] - W11
@fastmath @inbounds W11 += (W12 - W11)*(v.c[1]*v.corrRI - v.u[1])*dt*v.invDx # MODIFIED THIS LINE: added v.corrRI
@fastmath @inbounds W21 = 2.0*v.Q[1]/(v.A[1]-v.Ac) - W11 # IS THIS CORRECT? OR Q[1]/(v.A[1]-v.Ac)??

v.u[1], v.c[1] = inverseRiemannInvariants(W11, W21)
v.u[1], v.c[1] = inverseRiemannInvariants(W11, W21, v) # MODIFIED THIS LINE

if h.inlet_type == "Q"
@fastmath @inbounds v.A[1] = v.Q[1]/v.u[1]
@fastmath @inbounds v.A[1] = v.Q[1]/v.u[1] + v.Ac # MODIFIED THIS LINE: ADDED +v.Ac
v.P[1] = pressure(v.A[1], v.A0[1], v.beta[1], v.Pext)
else
v.A[1] = areaFromPressure(v.P[1], v.A0[1], v.beta[1], v.Pext)
@fastmath @inbounds v.Q[1] = v.u[1]*v.A[1]
@fastmath @inbounds v.Q[1] = v.u[1]*(v.A[1] - v.Ac) # MODIFIED THIS LINE: ADDED -v.Ac
end

end
Expand All @@ -93,8 +93,8 @@ end
Calculate Riemann invariants at the node `i` from `u` and `c`.
"""
function riemannInvariants(i :: Int, v :: Vessel)
@fastmath @inbounds W1 = v.u[i] - 4.0*v.c[i]
@fastmath @inbounds W2 = v.u[i] + 4.0*v.c[i]
@fastmath @inbounds W1 = v.u[i] - 4.0*v.c[i]*v.corrRI # MODIFIED THIS LINE
@fastmath @inbounds W2 = v.u[i] + 4.0*v.c[i]*v.corrRI # MODIFIED THIS LINE

return W1, W2
end
Expand All @@ -105,9 +105,9 @@ end

Calculate `u` and `c` given `W1` and `W2`
"""
function inverseRiemannInvariants(W1 :: Float64, W2 :: Float64)
function inverseRiemannInvariants(W1 :: Float64, W2 :: Float64, v :: Vessel)
@fastmath u = 0.5*(W1 + W2)
@fastmath c = (W2 - W1)*0.125
@fastmath c = (W2 - W1)*0.125/v.corrRI

return u, c
end
Expand Down Expand Up @@ -137,7 +137,7 @@ function setOutletBC(dt :: Float64, v :: Vessel)
outletCompatibility(dt, v)
elseif v.outlet == "wk3"
threeElementWindkessel(dt, v)
end
end
end


Expand All @@ -154,8 +154,8 @@ function outletCompatibility(dt :: Float64, v :: Vessel)
W2M += (W2M1 - W2M)*(v.u[end] + v.c[end])*dt/v.dx
W1M = v.W1M0 - v.Rt * (W2M - v.W2M0)

v.u[end], v.c[end] = inverseRiemannInvariants(W1M, W2M)
v.Q[end] = v.A[end]*v.u[end]
v.u[end], v.c[end] = inverseRiemannInvariants(W1M, W2M, v) # MODIFIED THIS LINE
v.Q[end] = (v.A[end] - v.Ac)*v.u[end] # MODIFIED THIS LINE: ADDED -v.Ac
end


Expand All @@ -169,10 +169,10 @@ solution of a Riemann problem at the 0D/1D interface.
function threeElementWindkessel(dt :: Float64, v :: Vessel)
Pout = 0.0

Al = v.A[end]
Al = v.A[end]
ul = v.u[end]

v.Pc += dt/v.Cc * (Al*ul - (v.Pc - Pout)/v.R2)
v.Pc += dt/v.Cc * (Al*ul - (v.Pc - Pout)/v.R2) # PROBABLY MUST BE (Al-v.Ac[end])*ul

As = Al

Expand All @@ -194,7 +194,7 @@ function threeElementWindkessel(dt :: Float64, v :: Vessel)
throw(e)
end

us = (pressure(As, v.A0[end], v.beta[end], v.Pext) - Pout)/(As*v.R1)
us = (pressure(As + v.Ac, v.A0[end], v.beta[end], v.Pext) - Pout)/(As*v.R1) # MODIFIED THIS LINE

v.A[end] = As
v.u[end] = us
Expand Down
38 changes: 24 additions & 14 deletions src/conjunctions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,14 +44,16 @@ end
Return the Jacobian for conjunction equations.
"""
function calculateJacobianConjunction(v1 :: Vessel, v2 :: Vessel, U, k)
@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[3]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[4]^4) # MODIFIED THIS LINE
@fastmath @inbounds U33 = U[3]*U[3]*U[3]
@fastmath @inbounds U43 = U[4]*U[4]*U[4]

@fastmath @inbounds J13 = 4.0*k[1]
@fastmath @inbounds J24 = -4.0*k[2]
@fastmath @inbounds J13 = 4.0*k[1]*(v1.Ac/(2*U[3]*g1) + g1) # MODIFIED THIS LINE
@fastmath @inbounds J24 = -4.0*k[2]*(v2.Ac/(2*U[4]*g2) + g2) # MODIFIED THIS LINE

@fastmath @inbounds J31 = U33*U[3]
@fastmath @inbounds J32 = -U43*U[4]
@fastmath @inbounds J31 = U33*U[3] - v1.Ac # MODIFIED THIS LINE
@fastmath @inbounds J32 = -U43*U[4] + v2.Ac # MODIFIED THIS LINE
@fastmath @inbounds J33 = 4.0*U[1]*U33
@fastmath @inbounds J34 = -4.0*U[2]*U43

Expand All @@ -72,9 +74,14 @@ end

Return the Riemann invariants at the conjunction node.
"""
function calculateWstarConjunction(U, k)
@fastmath @inbounds W1 = U[1] + 4.0*k[1]*U[3]
@fastmath @inbounds W2 = U[2] - 4.0*k[2]*U[4]
function calculateWstarConjunction(U, k, vessels :: Array{Vessel,1})
# MODIFIED THIS FUNCTION
v1 = vessels[1] # MODIFIED THIS LINE
v2 = vessels[2] # MODIFIED THIS LINE
@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[3]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[4]^4) # MODIFIED THIS LINE
@fastmath @inbounds W1 = U[1] + 4.0*k[1]*U[3]*g1 # MODIFIED THIS LINE
@fastmath @inbounds W2 = U[2] - 4.0*k[2]*U[4]*g2 # MODIFIED THIS LINE

return @SArray [W1, W2]
end
Expand All @@ -87,14 +94,17 @@ function calculateFConjunction(vessels :: Array{Vessel,1}, U, k, W)
v1 = vessels[1]
v2 = vessels[2]

@fastmath @inbounds g1 = sqrt(1-v1.Ac/U[3]^4) # MODIFIED THIS LINE
@fastmath @inbounds g2 = sqrt(1-v2.Ac/U[4]^4) # MODIFIED THIS LINE

@fastmath @inbounds U32 = U[3]*U[3]
@fastmath @inbounds U42 = U[4]*U[4]

@fastmath @inbounds f1 = U[1] + 4.0*k[1]*U[3] - W[1]
@fastmath @inbounds f1 = U[1] + 4.0*k[1]*U[3]*g1 - W[1] # MODIFIED THIS LINE

@fastmath @inbounds f2 = U[2] - 4.0*k[2]*U[4] - W[2]
@fastmath @inbounds f2 = U[2] - 4.0*k[2]*U[4]*g2 - W[2] # MODIFIED THIS LINE

@fastmath @inbounds f3 = U[1]*U32*U32 - U[2]*U42*U42
@fastmath @inbounds f3 = U[1]*(U32*U32 - v1.Ac) - U[2]*(U42*U42 - v2.Ac) # MODIFIED THIS LINE

f4 = 0.5*k[3]*U[1]*U[1] + v1.beta[end]*(U32*v1.s_inv_A0[end] - 1.0) -
(0.5*k[3]*U[2]*U[2] + v2.beta[1]*(U42*v2.s_inv_A0[1] - 1.0))
Expand All @@ -113,14 +123,14 @@ function updateConjunction(U, v1 :: Vessel, v2 :: Vessel)
@inbounds v2.u[1] = U[2]

@fastmath @inbounds v1.A[end] = U[3]*U[3]*U[3]*U[3]
@fastmath @inbounds v1.Q[end] = v1.u[end]*v1.A[end]
@fastmath @inbounds v1.Q[end] = v1.u[end]*(v1.A[end] - v1.Ac) # MODIFIED THIS LINE

@fastmath @inbounds v2.A[1] = U[4]*U[4]*U[4]*U[4]
@fastmath @inbounds v2.Q[1] = v2.u[1]*v2.A[1]
@fastmath @inbounds v2.Q[1] = v2.u[1]*(v2.A[1] - v2.Ac) # MODIFIED THIS LINE

@inbounds v1.P[end] = pressure(v1.A[end], v1.A0[end], v1.beta[end], v1.Pext)
@inbounds v2.P[1] = pressure(v2.A[1], v2.A0[1], v2.beta[1], v2.Pext)

@fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end])
@fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1])
@fastmath @inbounds v1.c[end] = waveSpeed(v1.A[end], v1.gamma[end], v1.Ac) # MODIFIED THIS LINE
@fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1], v2.Ac) # MODIFIED THIS LINE
end
Loading