From 98bbb6d6e4262fb764dd1d59efa23d6d47cce662 Mon Sep 17 00:00:00 2001 From: ibenemerito88 Date: Fri, 20 Dec 2019 14:48:26 +0000 Subject: [PATCH 1/4] first commit catheter branch --- .gitignore | 3 +- Project.toml | 1 + src/boundary_conditions.jl | 26 +-- src/conjunctions.jl | 36 ++-- src/initialise.jl | 35 +++- src/junctions.jl | 4 +- src/openBF.jl | 7 + src/solver.jl | 16 +- test/jointcat/.conjunction.yml.swp | Bin 0 -> 12288 bytes test/jointcat/.script.jl.swp | Bin 0 -> 12288 bytes test/jointcat/backupres/P_A.last | 100 ++++++++++ test/jointcat/backupres/P_A.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/P_P.last | 100 ++++++++++ test/jointcat/backupres/P_P.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/P_Q.last | 100 ++++++++++ test/jointcat/backupres/P_Q.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/P_c.last | 100 ++++++++++ test/jointcat/backupres/P_c.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/P_u.last | 100 ++++++++++ test/jointcat/backupres/P_u.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/d1_A.last | 100 ++++++++++ test/jointcat/backupres/d1_A.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/d1_P.last | 100 ++++++++++ test/jointcat/backupres/d1_P.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/d1_Q.last | 100 ++++++++++ test/jointcat/backupres/d1_Q.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/d1_c.last | 100 ++++++++++ test/jointcat/backupres/d1_c.out | 300 ++++++++++++++++++++++++++++ test/jointcat/backupres/d1_u.last | 100 ++++++++++ test/jointcat/backupres/d1_u.out | 300 ++++++++++++++++++++++++++++ test/jointcat/conjunction.yml | 44 ++++ test/jointcat/conjunction_inlet.dat | 100 ++++++++++ test/jointcat/plot.py | 18 ++ test/jointcat/script.jl | 21 ++ 34 files changed, 4267 insertions(+), 44 deletions(-) create mode 100644 test/jointcat/.conjunction.yml.swp create mode 100644 test/jointcat/.script.jl.swp create mode 100644 test/jointcat/backupres/P_A.last create mode 100644 test/jointcat/backupres/P_A.out create mode 100644 test/jointcat/backupres/P_P.last create mode 100644 test/jointcat/backupres/P_P.out create mode 100644 test/jointcat/backupres/P_Q.last create mode 100644 test/jointcat/backupres/P_Q.out create mode 100644 test/jointcat/backupres/P_c.last create mode 100644 test/jointcat/backupres/P_c.out create mode 100644 test/jointcat/backupres/P_u.last create mode 100644 test/jointcat/backupres/P_u.out create mode 100644 test/jointcat/backupres/d1_A.last create mode 100644 test/jointcat/backupres/d1_A.out create mode 100644 test/jointcat/backupres/d1_P.last create mode 100644 test/jointcat/backupres/d1_P.out create mode 100644 test/jointcat/backupres/d1_Q.last create mode 100644 test/jointcat/backupres/d1_Q.out create mode 100644 test/jointcat/backupres/d1_c.last create mode 100644 test/jointcat/backupres/d1_c.out create mode 100644 test/jointcat/backupres/d1_u.last create mode 100644 test/jointcat/backupres/d1_u.out create mode 100644 test/jointcat/conjunction.yml create mode 100644 test/jointcat/conjunction_inlet.dat create mode 100644 test/jointcat/plot.py create mode 100644 test/jointcat/script.jl diff --git a/.gitignore b/.gitignore index decf15d..18623ce 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ test/**/*_results/* **/.DS_Store cvs_results/* docs/build/* -Manifest.toml \ No newline at end of file +Manifest.toml +*.sh diff --git a/Project.toml b/Project.toml index 13e8f6b..225c99e 100644 --- a/Project.toml +++ b/Project.toml @@ -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" diff --git a/src/boundary_conditions.jl b/src/boundary_conditions.jl index a69e343..8721c18 100644 --- a/src/boundary_conditions.jl +++ b/src/boundary_conditions.jl @@ -72,16 +72,16 @@ function inletCompatibility(dt :: Float64, v :: Vessel, h :: Heart) 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 W21 = 2.0*v.Q[1]/(v.A[1]-v.Ac[1]) - W11 # IS THIS CORRECT? OR Q[1]/(v.A[1]-v.Ac[1])?? - 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[1] # MODIFIED THIS LINE: ADDED +v.Ac[1] 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[1]) # MODIFIED THIS LINE: ADDED -v.Ac[1] end end @@ -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 @@ -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 @@ -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[end])*v.u[end] # MODIFIED THIS LINE: ADDED -v.Ac[end] end @@ -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] - v.Ac[end] # MODIFIED THIS LINE 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 @@ -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[end], v.A0[end], v.beta[end], v.Pext) - Pout)/(As*v.R1) # MODIFIED THIS LINE v.A[end] = As v.u[end] = us diff --git a/src/conjunctions.jl b/src/conjunctions.jl index fab17bd..a410444 100644 --- a/src/conjunctions.jl +++ b/src/conjunctions.jl @@ -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[end]/U[3]^4) # MODIFIED THIS LINE + @fastmath @inbounds g2 = sqrt(1-v2.Ac[1]/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[end]/(2*U[3]*g1) + g1) # MODIFIED THIS LINE + @fastmath @inbounds J24 = -4.0*k[2]*(v2.Ac[1]/(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[end] # MODIFIED THIS LINE + @fastmath @inbounds J32 = -U43*U[4] + v2.Ac[1] # MODIFIED THIS LINE @fastmath @inbounds J33 = 4.0*U[1]*U33 @fastmath @inbounds J34 = -4.0*U[2]*U43 @@ -72,9 +74,12 @@ 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, v1:: Vessel, v2 :: Vessel) + # MODIFIED THIS FUNCTION + @fastmath @inbounds g1 = sqrt(1-v1.Ac[end]/U[3]^4) # MODIFIED THIS LINE + @fastmath @inbounds g2 = sqrt(1-v2.Ac[1]/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 @@ -87,14 +92,17 @@ function calculateFConjunction(vessels :: Array{Vessel,1}, U, k, W) v1 = vessels[1] v2 = vessels[2] + @fastmath @inbounds g1 = sqrt(1-v1.Ac[end]/U[3]^4) # MODIFIED THIS LINE + @fastmath @inbounds g2 = sqrt(1-v2.Ac[1]/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[end]) - U[2]*(U42*U42 - v2.Ac[1]) # 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)) @@ -113,14 +121,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[end]) # 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[1]) # 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[end]) # MODIFIED THIS LINE + @fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1], v2.Ac[1]) # MODIFIED THIS LINE end diff --git a/src/initialise.jl b/src/initialise.jl index 417b884..a680712 100644 --- a/src/initialise.jl +++ b/src/initialise.jl @@ -244,7 +244,7 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju L = vessel_data["L"] E = vessel_data["E"] - Rp, Rd = computeRadii(vessel_data) + Rp, Rd, Rc = computeRadii(vessel_data) # MODIFIED THIS LINE Pext = getPext(vessel_data) M, dx, invDx, halfDx, invDxSq = meshVessel(vessel_data, L) h0 = initialiseThickness(vessel_data, M) @@ -256,6 +256,7 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju Q = zeros(Float64, M) P = zeros(Float64, M) A = zeros(Float64, M) + Ac = zeros(Float64, M) # MODIFIED THIS LINE u = zeros(Float64, M) c = zeros(Float64, M) A0 = zeros(Float64, M) @@ -318,12 +319,13 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju inv_A0[i] = 1.0/A0[i] s_inv_A0[i] = sqrt(inv_A0[i]) A[i] = A0[i] + Ac[i] = pi*Rc*Rc # MODIFIED THIS LINE beta[i] = s_inv_A0[i]*h0*s_pi_E_over_sigma_squared gamma[i] = beta[i]*one_over_rho_s_p/R0[i] s_15_gamma[i] = sqrt(1.5*gamma[i]) gamma_ghost[i+1] = gamma[i] P[i] = pressure(1.0, beta[i], Pext) - c[i] = waveSpeed(A[i], gamma[i]) + c[i] = waveSpeed(A[i], gamma[i], Ac[i]) wallE[i] = 3.0*beta[i]*radius_slope*inv_A0[i]*s_pi*blood.rho_inv if phi != 0.0 wallVb[i] = Cv*s_inv_A0[i]*invDxSq @@ -349,8 +351,15 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju UM1Q = 0.0 UM2Q = 0.0 - W1M0 = u[end] - 4.0*c[end] - W2M0 = u[end] + 4.0*c[end] + # Compute correction factor to Riemann Invariants + if haskey(vessel_data, "Rc") + corrRI = (A[end]/Ac[end])^0.25*drummond2F1(0.25,0.5,1.5,1-A[end]/Ac[end])/2 + else + corrRI = 1 + end + + W1M0 = u[end] - 4.0*c[end]*corrRI # MODIFIED THIS LINE + W2M0 = u[end] + 4.0*c[end]*corrRI # MODIFIED THIS LINE node2 = convert(Int, floor(M*0.25)) node3 = convert(Int, floor(M*0.5)) @@ -390,6 +399,7 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju slope, flux, uStar, vA, vQ, dU, slopesA, slopesQ, Al, Ar, Ql, Qr, Fl, Fr, + Ac, corrRI, outlet) end @@ -438,13 +448,18 @@ If only a constant lumen radius is defined, return the same value for proximal a distal variables, `Rp` and `Rd`, respectively. """ function computeRadii(vessel :: Dict{Any,Any}) + if ~haskey(vessel, "Rc") # Modified from there: first define Rc + Rc = 0 + else + Rc = vessel["Rc"] + end # to here if ~haskey(vessel, "R0") Rp = vessel["Rp"] Rd = vessel["Rd"] - return Rp, Rd + return Rp, Rd, Rc # MODIFIED THIS LINE else R0 = vessel["R0"] - return R0, R0 + return R0, R0, Rc # MODIFIED THIS LINE end end @@ -570,12 +585,16 @@ where gamma_v (`gamma_profile`) is either specified in the vessel definition or assumed equal to `9` (plug-flow). """ function computeViscousTerm(vessel_data :: Dict{Any,Any}, blood :: Blood) - if haskey(vessel_data, "gamma_profile") + + if haskey(vessel_data,"Rc") && (vessel_data["Rc"])>0 # MODIFIED THIS LINE + return 2*blood.mu*blood.rho_inv*sqrt(15) + elseif haskey(vessel_data, "gamma_profile") gamma_profile = vessel_data["gamma_profile"] else gamma_profile = 9 - end return 2*(gamma_profile + 2)*pi*blood.mu*blood.rho_inv + end + end diff --git a/src/junctions.jl b/src/junctions.jl index 03b86c2..dcd648e 100644 --- a/src/junctions.jl +++ b/src/junctions.jl @@ -37,7 +37,7 @@ end Solve the non-linear junction system by means of Newton-Raphson method. """ function newtonRaphson(vessels :: Array{Vessel,1}, J, U, k, funW, funF) - W = funW(U, k) + W = funW(U, k, vessels[1], vessels[2]) # MODIFIED THIS LINE F = funF(vessels, U, k, W) nr_toll_U = 1e-5 @@ -71,7 +71,7 @@ function newtonRaphson(vessels :: Array{Vessel,1}, J, U, k, funW, funF) else U = U_new - W = funW(U, k) + W = funW(U, k, vessels[1], vessels[2]) # MODIFIED THIS LINE F = funF(vessels, U, k, W) end end diff --git a/src/openBF.jl b/src/openBF.jl index 31d5990..f63895d 100644 --- a/src/openBF.jl +++ b/src/openBF.jl @@ -25,6 +25,9 @@ module openBF using DelimitedFiles using LinearAlgebra using Printf + # IMPORT PACKAGE FOR HYPERGEOMETRIC FUNCTIONS + using HypergeometricFunctions + import HypergeometricFunctions: drummond2F1 """ @@ -167,6 +170,10 @@ module openBF Fl :: Array{Float64,2} Fr :: Array{Float64,2} + # Catheter + Ac :: Array{Float64,1} # Area of the catheter + corrRI :: Float64 # Correction factor to Riemann Invariants + #Outlet type outlet :: String end diff --git a/src/solver.jl b/src/solver.jl index a91e913..6b9f843 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -36,8 +36,8 @@ end Compute pulse wave velocity at the given node. """ -function waveSpeed(A :: Float64, gamma :: Float64) - return sqrt(3*gamma*sqrt(A)*0.5) +function waveSpeed(A :: Float64, gamma :: Float64, Ac :: Float64) + return sqrt(3*gamma*0.5*sqrt(A))*sqrt(1-Ac/A) # MODIFIED THIS LINE end function waveSpeedSA(sA :: Float64, gamma :: Float64) @@ -216,11 +216,15 @@ function muscl(v :: Vessel, dt :: Float64, b :: Blood) #source term @fastmath @inbounds @simd for i = 1:v.M s_A = sqrt(v.A[i]) - Si = - v.viscT*v.Q[i]/v.A[i] - v.wallE[i]*(s_A - v.s_A0[i])*v.A[i] + if v.Ac[1]>0 + Si = -v.viscT*v.Q[i]/(s_A-sqrt(v.Ac[i]))^2; + else + Si = - v.viscT*v.Q[i]/v.A[i] - v.wallE[i]*(s_A - v.s_A0[i])*v.A[i] + end v.Q[i] += dt*Si v.P[i] = pressure(s_A*v.s_inv_A0[i], v.beta[i], v.Pext) - v.c[i] = waveSpeedSA(s_A, v.gamma[i]) + v.c[i] = waveSpeed(v.A[i], v.gamma[i], v.Ac[i]) end @@ -241,7 +245,7 @@ function muscl(v :: Vessel, dt :: Float64, b :: Blood) end @fastmath @inbounds @simd for i = 1:v.M - v.u[i] = v.Q[i]/v.A[i] + v.u[i] = v.Q[i]/(v.A[i] - v.Ac[i]) # MODIFIED THIS LINE end end @@ -254,7 +258,7 @@ function computeFlux(v :: Vessel, A :: Array{Float64,1}, Q :: Array{Float64,1}, Flux :: Array{Float64,2}) @fastmath @inbounds @simd for i in 1:v.M+2 Flux[1,i] = Q[i] - Flux[2,i] = Q[i]*Q[i]/A[i] + v.gamma_ghost[i]*A[i]*sqrt(A[i]) + Flux[2,i] = Q[i]*Q[i]/(A[i] - v.Ac[i]) + v.gamma_ghost[i]*A[i]*sqrt(A[i]) - 3*v.gamma_ghost[i]*v.Ac[i]*sqrt(A[i]) end return Flux diff --git a/test/jointcat/.conjunction.yml.swp b/test/jointcat/.conjunction.yml.swp new file mode 100644 index 0000000000000000000000000000000000000000..0a28cc2e0088200a701d5907365c5ff605b1fd93 GIT binary patch literal 12288 zcmeI2&x;gC6vt~Nxu~0;dVe|xS9E%MW>+2i7FPwqz$!!tVO?svYi4`4x@zdInsq_V zNkKg6(W8Gr4!Ic!5Rd)^IpkYCy|@XeIe1C(3O;mq^{e-)K3#K~8pKyF zd`bN%7HFRcaeVua)$ixN6MK&dkx!)w-S284uI&fjIRDnrn<21swzQLhz1CT=mCJ0} zS(as%+o`NluVvY>6J^ev@7Ocv*JM`YW~alc%5}yzH!hW~lUg(OCUv|~T9#k3 z0XE>>u-%^*i@k3BQ*RwUM4x~56?<~18m@bWWY`4#J5=Qo@UL* zo9o=)uFWG~*Z><~18jf|umLu}2G{@_U;}J`4X}ZKp#hl*@oJwCKkUck!~g%y8eV=Z z#P47m{0Q!Yd*BWjgD=3R;Q2>FJOjUhU%@ZnA-DmqgEep!#Nb2l0oa}w;t6;RegY4` zZEy>0f+^5o03~3Dk^N|wXbMl zMIN0EyNqF{qy8xjL-H0)N_#CYq$cePWvG-{9^#88{9;1gUesMu?H+=)G9$PC`mHkw z&8k

xBw!(sH15LRoiKFgKV!8D9)kiDdO$XMMw5tAutPd0wmgqfENjQ3+jo<9e6= z$$C>vEL!T%){7;F)*>9XURO?)&4NYfeeY^|AuEbydF5bu!O3!@+)*+%QzESaqza*wZlyY8p@E`GPE{82hrVhyqD6(oHt|LCA;}J#_Xu20hxp!c zbh=KwjM+=G^TkDFHD*~)xbO)rMLQcW+cL7qQ<~;gYOCCB(Kjpd#ob%gDpHo~HWWUI i{T%!$GYi-lb&9Ktt74$7&5#qDcZ}ZzA=6)XlK2zOUp%q^ literal 0 HcmV?d00001 diff --git a/test/jointcat/.script.jl.swp b/test/jointcat/.script.jl.swp new file mode 100644 index 0000000000000000000000000000000000000000..650c0db71384387b151ef75b70edcafe28e52554 GIT binary patch literal 12288 zcmeI2v2GJV5QZnz1P~GZaY)Dp$LD~AXk3s)jwk{s7z(VDwZ0wfOZIM0yL*Ndz#Gu< z8mW+Y1|)jkfQAwxQ8Ifzlf*%x4w{pGuXc8JcIMl2VR=TK=MT4o`|S?Ua|Ph>;jj1; z9Khx}Kspts#*_Z5myS^y}z}wod9K+7L#jdMa#EDq$jBV7}g) z4-T!15}ogsdE!!?M>a8O;o7_E%#(PL2{3`v1eW3Mowe$#ej~gdtXzF|+J!AlfC(@G zCcp%k025#WOn?dew**|V1n;Q23w4Fp>TzjdtfxFM0Vco%m;e)C0!)AjFaajO1egF5 zU;^imfRF%RF9K{WQ}Xox|NQs=kIMkxNFPaWNPW^h()UXMUq}a}52PXKCTWGVL^`B6 zKS`fS?@1Zy1#{If|1f%wkTLD zIDSAAqzk%1eLy0VZ>Dy*iOPxoQIL(~rD_5sntHX-(1tZA1SKtp$$p~HlBIJpaWYVE zxUZqmY3{s-5&0M!omOC-F1iyG# Date: Thu, 26 Dec 2019 21:07:25 +0000 Subject: [PATCH 2/4] Ac scalar --- .gitignore | 1 + src/boundary_conditions.jl | 14 +++++++------- src/conjunctions.jl | 30 ++++++++++++++--------------- src/initialise.jl | 19 +++++++++--------- src/openBF.jl | 6 +++--- src/solver.jl | 16 ++++++--------- test/jointcat/.conjunction.yml.swp | Bin 12288 -> 0 bytes test/jointcat/.script.jl.swp | Bin 12288 -> 0 bytes test/jointcat/conjunction.yml | 4 ++-- 9 files changed, 43 insertions(+), 47 deletions(-) delete mode 100644 test/jointcat/.conjunction.yml.swp delete mode 100644 test/jointcat/.script.jl.swp diff --git a/.gitignore b/.gitignore index 18623ce..93153e6 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ cvs_results/* docs/build/* Manifest.toml *.sh +*.swp diff --git a/src/boundary_conditions.jl b/src/boundary_conditions.jl index 8721c18..683ba94 100644 --- a/src/boundary_conditions.jl +++ b/src/boundary_conditions.jl @@ -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]-v.Ac[1]) - W11 # IS THIS CORRECT? OR Q[1]/(v.A[1]-v.Ac[1])?? + @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) # MODIFIED THIS LINE if h.inlet_type == "Q" - @fastmath @inbounds v.A[1] = v.Q[1]/v.u[1] + v.Ac[1] # MODIFIED THIS LINE: ADDED +v.Ac[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] - v.Ac[1]) # MODIFIED THIS LINE: ADDED -v.Ac[1] + @fastmath @inbounds v.Q[1] = v.u[1]*(v.A[1] - v.Ac) # MODIFIED THIS LINE: ADDED -v.Ac end end @@ -155,7 +155,7 @@ function outletCompatibility(dt :: Float64, v :: Vessel) W1M = v.W1M0 - v.Rt * (W2M - v.W2M0) v.u[end], v.c[end] = inverseRiemannInvariants(W1M, W2M, v) # MODIFIED THIS LINE - v.Q[end] = (v.A[end] - v.Ac[end])*v.u[end] # MODIFIED THIS LINE: ADDED -v.Ac[end] + v.Q[end] = (v.A[end] - v.Ac)*v.u[end] # MODIFIED THIS LINE: ADDED -v.Ac end @@ -169,7 +169,7 @@ solution of a Riemann problem at the 0D/1D interface. function threeElementWindkessel(dt :: Float64, v :: Vessel) Pout = 0.0 - Al = v.A[end] - v.Ac[end] # MODIFIED THIS LINE + Al = v.A[end] ul = v.u[end] v.Pc += dt/v.Cc * (Al*ul - (v.Pc - Pout)/v.R2) # PROBABLY MUST BE (Al-v.Ac[end])*ul @@ -194,7 +194,7 @@ function threeElementWindkessel(dt :: Float64, v :: Vessel) throw(e) end - us = (pressure(As + v.Ac[end], v.A0[end], v.beta[end], v.Pext) - Pout)/(As*v.R1) # MODIFIED THIS LINE + 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 diff --git a/src/conjunctions.jl b/src/conjunctions.jl index a410444..786d883 100644 --- a/src/conjunctions.jl +++ b/src/conjunctions.jl @@ -44,16 +44,16 @@ end Return the Jacobian for conjunction equations. """ function calculateJacobianConjunction(v1 :: Vessel, v2 :: Vessel, U, k) - @fastmath @inbounds g1 = sqrt(1-v1.Ac[end]/U[3]^4) # MODIFIED THIS LINE - @fastmath @inbounds g2 = sqrt(1-v2.Ac[1]/U[4]^4) # 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 U33 = U[3]*U[3]*U[3] @fastmath @inbounds U43 = U[4]*U[4]*U[4] - @fastmath @inbounds J13 = 4.0*k[1]*(v1.Ac[end]/(2*U[3]*g1) + g1) # MODIFIED THIS LINE - @fastmath @inbounds J24 = -4.0*k[2]*(v2.Ac[1]/(2*U[4]*g2) + g2) # MODIFIED THIS LINE + @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] - v1.Ac[end] # MODIFIED THIS LINE - @fastmath @inbounds J32 = -U43*U[4] + v2.Ac[1] # MODIFIED THIS LINE + @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 @@ -76,8 +76,8 @@ Return the Riemann invariants at the conjunction node. """ function calculateWstarConjunction(U, k, v1:: Vessel, v2 :: Vessel) # MODIFIED THIS FUNCTION - @fastmath @inbounds g1 = sqrt(1-v1.Ac[end]/U[3]^4) # MODIFIED THIS LINE - @fastmath @inbounds g2 = sqrt(1-v2.Ac[1]/U[4]^4) # 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 @@ -92,8 +92,8 @@ function calculateFConjunction(vessels :: Array{Vessel,1}, U, k, W) v1 = vessels[1] v2 = vessels[2] - @fastmath @inbounds g1 = sqrt(1-v1.Ac[end]/U[3]^4) # MODIFIED THIS LINE - @fastmath @inbounds g2 = sqrt(1-v2.Ac[1]/U[4]^4) # 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 U32 = U[3]*U[3] @fastmath @inbounds U42 = U[4]*U[4] @@ -102,7 +102,7 @@ function calculateFConjunction(vessels :: Array{Vessel,1}, U, k, W) @fastmath @inbounds f2 = U[2] - 4.0*k[2]*U[4]*g2 - W[2] # MODIFIED THIS LINE - @fastmath @inbounds f3 = U[1]*(U32*U32 - v1.Ac[end]) - U[2]*(U42*U42 - v2.Ac[1]) # MODIFIED THIS LINE + @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)) @@ -121,14 +121,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] - v1.Ac[end]) # MODIFIED THIS LINE + @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] - v2.Ac[1]) # MODIFIED THIS LINE + @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], v1.Ac[end]) # MODIFIED THIS LINE - @fastmath @inbounds v2.c[1] = waveSpeed(v2.A[1], v2.gamma[1], v2.Ac[1]) # MODIFIED THIS LINE + @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 diff --git a/src/initialise.jl b/src/initialise.jl index a680712..5b61f5d 100644 --- a/src/initialise.jl +++ b/src/initialise.jl @@ -256,7 +256,7 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju Q = zeros(Float64, M) P = zeros(Float64, M) A = zeros(Float64, M) - Ac = zeros(Float64, M) # MODIFIED THIS LINE + Ac = pi*Rc*Rc # MODIFIED THIS LINE u = zeros(Float64, M) c = zeros(Float64, M) A0 = zeros(Float64, M) @@ -319,13 +319,12 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju inv_A0[i] = 1.0/A0[i] s_inv_A0[i] = sqrt(inv_A0[i]) A[i] = A0[i] - Ac[i] = pi*Rc*Rc # MODIFIED THIS LINE beta[i] = s_inv_A0[i]*h0*s_pi_E_over_sigma_squared gamma[i] = beta[i]*one_over_rho_s_p/R0[i] s_15_gamma[i] = sqrt(1.5*gamma[i]) gamma_ghost[i+1] = gamma[i] P[i] = pressure(1.0, beta[i], Pext) - c[i] = waveSpeed(A[i], gamma[i], Ac[i]) + c[i] = waveSpeed(A[i], gamma[i], Ac) wallE[i] = 3.0*beta[i]*radius_slope*inv_A0[i]*s_pi*blood.rho_inv if phi != 0.0 wallVb[i] = Cv*s_inv_A0[i]*invDxSq @@ -353,9 +352,9 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju # Compute correction factor to Riemann Invariants if haskey(vessel_data, "Rc") - corrRI = (A[end]/Ac[end])^0.25*drummond2F1(0.25,0.5,1.5,1-A[end]/Ac[end])/2 + corrRI = (A[end]/Ac)^0.25*drummond2F1(0.25,0.5,1.5,1-A[end]/Ac)/2 else - corrRI = 1 + corrRI = 1.0 end W1M0 = u[end] - 4.0*c[end]*corrRI # MODIFIED THIS LINE @@ -399,7 +398,7 @@ function buildVessel(ID :: Int, vessel_data :: Dict{Any,Any}, blood :: Blood, ju slope, flux, uStar, vA, vQ, dU, slopesA, slopesQ, Al, Ar, Ql, Qr, Fl, Fr, - Ac, corrRI, + Ac, corrRI, outlet) end @@ -449,10 +448,10 @@ distal variables, `Rp` and `Rd`, respectively. """ function computeRadii(vessel :: Dict{Any,Any}) if ~haskey(vessel, "Rc") # Modified from there: first define Rc - Rc = 0 + Rc = 0 else - Rc = vessel["Rc"] - end # to here + Rc = vessel["Rc"] + end # to here if ~haskey(vessel, "R0") Rp = vessel["Rp"] Rd = vessel["Rd"] @@ -586,7 +585,7 @@ assumed equal to `9` (plug-flow). """ function computeViscousTerm(vessel_data :: Dict{Any,Any}, blood :: Blood) - if haskey(vessel_data,"Rc") && (vessel_data["Rc"])>0 # MODIFIED THIS LINE + if haskey(vessel_data,"Rc") # MODIFIED THIS LINE return 2*blood.mu*blood.rho_inv*sqrt(15) elseif haskey(vessel_data, "gamma_profile") gamma_profile = vessel_data["gamma_profile"] diff --git a/src/openBF.jl b/src/openBF.jl index f63895d..307eb68 100644 --- a/src/openBF.jl +++ b/src/openBF.jl @@ -170,9 +170,9 @@ module openBF Fl :: Array{Float64,2} Fr :: Array{Float64,2} - # Catheter - Ac :: Array{Float64,1} # Area of the catheter - corrRI :: Float64 # Correction factor to Riemann Invariants + # Catheter + Ac :: Float64 # Area of the catheter + corrRI :: Float64 # Correction factor to Riemann Invariants #Outlet type outlet :: String diff --git a/src/solver.jl b/src/solver.jl index 6b9f843..028771f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -32,9 +32,9 @@ end """ - waveSpeed(A :: Float64, gamma :: Float64) + waveSpeed(A :: Float64, gamma :: Float64, Ac :: Float64) -Compute pulse wave velocity at the given node. +Compute pulse wave velocity at the given node. It includes the catheter. """ function waveSpeed(A :: Float64, gamma :: Float64, Ac :: Float64) return sqrt(3*gamma*0.5*sqrt(A))*sqrt(1-Ac/A) # MODIFIED THIS LINE @@ -216,15 +216,11 @@ function muscl(v :: Vessel, dt :: Float64, b :: Blood) #source term @fastmath @inbounds @simd for i = 1:v.M s_A = sqrt(v.A[i]) - if v.Ac[1]>0 - Si = -v.viscT*v.Q[i]/(s_A-sqrt(v.Ac[i]))^2; - else - Si = - v.viscT*v.Q[i]/v.A[i] - v.wallE[i]*(s_A - v.s_A0[i])*v.A[i] - end + Si = - v.viscT*v.Q[i]/(s_A-sqrt(v.Ac))^2 - v.wallE[i]*(s_A - v.s_A0[i])*v.A[i] # MODIFIED THIS LINE v.Q[i] += dt*Si v.P[i] = pressure(s_A*v.s_inv_A0[i], v.beta[i], v.Pext) - v.c[i] = waveSpeed(v.A[i], v.gamma[i], v.Ac[i]) + v.c[i] = waveSpeed(v.A[i], v.gamma[i], v.Ac) end @@ -245,7 +241,7 @@ function muscl(v :: Vessel, dt :: Float64, b :: Blood) end @fastmath @inbounds @simd for i = 1:v.M - v.u[i] = v.Q[i]/(v.A[i] - v.Ac[i]) # MODIFIED THIS LINE + v.u[i] = v.Q[i]/(v.A[i] - v.Ac) # MODIFIED THIS LINE end end @@ -258,7 +254,7 @@ function computeFlux(v :: Vessel, A :: Array{Float64,1}, Q :: Array{Float64,1}, Flux :: Array{Float64,2}) @fastmath @inbounds @simd for i in 1:v.M+2 Flux[1,i] = Q[i] - Flux[2,i] = Q[i]*Q[i]/(A[i] - v.Ac[i]) + v.gamma_ghost[i]*A[i]*sqrt(A[i]) - 3*v.gamma_ghost[i]*v.Ac[i]*sqrt(A[i]) + Flux[2,i] = Q[i]*Q[i]/(A[i] - v.Ac) + v.gamma_ghost[i]*A[i]*sqrt(A[i]) - 3*v.gamma_ghost[i]*v.Ac*sqrt(A[i]) # MODIFIED end return Flux diff --git a/test/jointcat/.conjunction.yml.swp b/test/jointcat/.conjunction.yml.swp deleted file mode 100644 index 0a28cc2e0088200a701d5907365c5ff605b1fd93..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12288 zcmeI2&x;gC6vt~Nxu~0;dVe|xS9E%MW>+2i7FPwqz$!!tVO?svYi4`4x@zdInsq_V zNkKg6(W8Gr4!Ic!5Rd)^IpkYCy|@XeIe1C(3O;mq^{e-)K3#K~8pKyF zd`bN%7HFRcaeVua)$ixN6MK&dkx!)w-S284uI&fjIRDnrn<21swzQLhz1CT=mCJ0} zS(as%+o`NluVvY>6J^ev@7Ocv*JM`YW~alc%5}yzH!hW~lUg(OCUv|~T9#k3 z0XE>>u-%^*i@k3BQ*RwUM4x~56?<~18m@bWWY`4#J5=Qo@UL* zo9o=)uFWG~*Z><~18jf|umLu}2G{@_U;}J`4X}ZKp#hl*@oJwCKkUck!~g%y8eV=Z z#P47m{0Q!Yd*BWjgD=3R;Q2>FJOjUhU%@ZnA-DmqgEep!#Nb2l0oa}w;t6;RegY4` zZEy>0f+^5o03~3Dk^N|wXbMl zMIN0EyNqF{qy8xjL-H0)N_#CYq$cePWvG-{9^#88{9;1gUesMu?H+=)G9$PC`mHkw z&8k

xBw!(sH15LRoiKFgKV!8D9)kiDdO$XMMw5tAutPd0wmgqfENjQ3+jo<9e6= z$$C>vEL!T%){7;F)*>9XURO?)&4NYfeeY^|AuEbydF5bu!O3!@+)*+%QzESaqza*wZlyY8p@E`GPE{82hrVhyqD6(oHt|LCA;}J#_Xu20hxp!c zbh=KwjM+=G^TkDFHD*~)xbO)rMLQcW+cL7qQ<~;gYOCCB(Kjpd#ob%gDpHo~HWWUI i{T%!$GYi-lb&9Ktt74$7&5#qDcZ}ZzA=6)XlK2zOUp%q^ diff --git a/test/jointcat/.script.jl.swp b/test/jointcat/.script.jl.swp deleted file mode 100644 index 650c0db71384387b151ef75b70edcafe28e52554..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12288 zcmeI2v2GJV5QZnz1P~GZaY)Dp$LD~AXk3s)jwk{s7z(VDwZ0wfOZIM0yL*Ndz#Gu< z8mW+Y1|)jkfQAwxQ8Ifzlf*%x4w{pGuXc8JcIMl2VR=TK=MT4o`|S?Ua|Ph>;jj1; z9Khx}Kspts#*_Z5myS^y}z}wod9K+7L#jdMa#EDq$jBV7}g) z4-T!15}ogsdE!!?M>a8O;o7_E%#(PL2{3`v1eW3Mowe$#ej~gdtXzF|+J!AlfC(@G zCcp%k025#WOn?dew**|V1n;Q23w4Fp>TzjdtfxFM0Vco%m;e)C0!)AjFaajO1egF5 zU;^imfRF%RF9K{WQ}Xox|NQs=kIMkxNFPaWNPW^h()UXMUq}a}52PXKCTWGVL^`B6 zKS`fS?@1Zy1#{If|1f%wkTLD zIDSAAqzk%1eLy0VZ>Dy*iOPxoQIL(~rD_5sntHX-(1tZA1SKtp$$p~HlBIJpaWYVE zxUZqmY3{s-5&0M!omOC-F1iyG# Date: Sat, 4 Jan 2020 01:18:20 +0000 Subject: [PATCH 3/4] add bifurcations, conjunctions --- src/anastomosis.jl | 56 ++++++---- src/bifurcations.jl | 54 ++++++---- src/conjunctions.jl | 10 +- src/junctions.jl | 4 +- test/jointbif/bifurcation.yml | 61 +++++++++++ test/jointbif/bifurcation_inlet.dat | 100 ++++++++++++++++++ test/jointbif/plot.py | 154 ++++++++++++++++++++++++++++ test/jointcat/conjunction.yml | 6 +- test/jointcat/script.jl | 4 +- 9 files changed, 396 insertions(+), 53 deletions(-) create mode 100644 test/jointbif/bifurcation.yml create mode 100644 test/jointbif/bifurcation_inlet.dat create mode 100644 test/jointbif/plot.py diff --git a/src/anastomosis.jl b/src/anastomosis.jl index 4fe435e..4191df4 100644 --- a/src/anastomosis.jl +++ b/src/anastomosis.jl @@ -46,19 +46,20 @@ end Return the Jacobian for anastomosis equations. """ function calculateJacobianAnastomosis(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 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 + @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 @@ -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 + @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 @@ -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)) @@ -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 diff --git a/src/bifurcations.jl b/src/bifurcations.jl index 199944d..bc31523 100644 --- a/src/bifurcations.jl +++ b/src/bifurcations.jl @@ -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 @@ -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 @@ -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)) @@ -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 diff --git a/src/conjunctions.jl b/src/conjunctions.jl index 786d883..d247d98 100644 --- a/src/conjunctions.jl +++ b/src/conjunctions.jl @@ -74,11 +74,13 @@ end Return the Riemann invariants at the conjunction node. """ -function calculateWstarConjunction(U, k, v1:: Vessel, v2 :: Vessel) +function calculateWstarConjunction(U, k, vessels :: Array{Vessel,1}) # MODIFIED THIS FUNCTION - @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 + 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] diff --git a/src/junctions.jl b/src/junctions.jl index dcd648e..1184552 100644 --- a/src/junctions.jl +++ b/src/junctions.jl @@ -37,7 +37,7 @@ end Solve the non-linear junction system by means of Newton-Raphson method. """ function newtonRaphson(vessels :: Array{Vessel,1}, J, U, k, funW, funF) - W = funW(U, k, vessels[1], vessels[2]) # MODIFIED THIS LINE + W = funW(U, k, vessels) # MODIFIED THIS LINE F = funF(vessels, U, k, W) nr_toll_U = 1e-5 @@ -71,7 +71,7 @@ function newtonRaphson(vessels :: Array{Vessel,1}, J, U, k, funW, funF) else U = U_new - W = funW(U, k, vessels[1], vessels[2]) # MODIFIED THIS LINE + W = funW(U, k, vessels) # MODIFIED THIS LINE F = funF(vessels, U, k, W) end end diff --git a/test/jointbif/bifurcation.yml b/test/jointbif/bifurcation.yml new file mode 100644 index 0000000..172e0bc --- /dev/null +++ b/test/jointbif/bifurcation.yml @@ -0,0 +1,61 @@ +project name: bifurcation + +blood: + rho: 1060.0 # density [kg/m^3] + mu: 4.e-3 # dynamic viscosity [Pa⋅s] + +solver: + Ccfl: 0.9 # Courant number + cycles: 100 # maximum number of cycles + jump: 100 # timesteps per cycle to be saved + convergence tolerance: 5.0 # percentage value + +network: + - label: P + # Int. carotid II (Melis et al., 2019) + sn: 1 # proximal (source) node + tn: 2 # distal (target) node + + L: 0.5e-2 # length (m) + + R0: 0.2e-2 # proximal lumen radius (m) + Rc: 0.15e-2 # catheter radius + + E: 1600.0e3 # Young's modulus (Pa) + h0: 0.82e-3 + + inlet: Q + inlet file: bifurcation_inlet.dat + inlet number: 1 + + - label: d1 + # ACA I + sn: 2 + tn: 3 + + L: 1.2e-2 # length (m) + + R0: 0.117e-2 # proximal lumen radius (m) + + E: 1600.0e3 # Young's modulus (Pa) + + outlet: wk3 + R1: 6.8123e7 + R2: 3.1013e9 + Cc: 3.6664e-10 + + - label: d2 + # MCA + sn: 2 + tn: 4 + + L: 11.9e-2 # length (m) + + R0: 0.143e-2 # proximal lumen radius (m) + + E: 1600.0e3 # Young's modulus (Pa) + + outlet: wk3 + R1: 6.8123e7 + R2: 3.1013e9 + Cc: 3.6664e-10 diff --git a/test/jointbif/bifurcation_inlet.dat b/test/jointbif/bifurcation_inlet.dat new file mode 100644 index 0000000..92ba01c --- /dev/null +++ b/test/jointbif/bifurcation_inlet.dat @@ -0,0 +1,100 @@ +0.000000000000000000e+00 -5.239489231023915536e-07 +1.111111111111111154e-02 -1.810043543947649782e-06 +2.222222222222222307e-02 -4.106752263203652509e-06 +3.333333333333333287e-02 -6.307155406759066627e-06 +4.444444444444444614e-02 -7.453316028919026914e-06 +5.555555555555555941e-02 -7.309722996010778137e-06 +6.666666666666666574e-02 -6.197508131371585217e-06 +7.777777777777777901e-02 -4.412937386033050817e-06 +8.888888888888889228e-02 -1.865784178552097541e-06 +1.000000000000000056e-01 1.765273247947420437e-06 +1.111111111111111188e-01 6.625078440213615815e-06 +1.222222222222222321e-01 1.252104753433897402e-05 +1.333333333333333315e-01 1.920573178717305051e-05 +1.444444444444444586e-01 2.677024308819679155e-05 +1.555555555555555580e-01 3.562926121035754368e-05 +1.666666666666666852e-01 4.598767697454575590e-05 +1.777777777777777846e-01 5.727420543835721757e-05 +1.888888888888888840e-01 6.813722078324621991e-05 +2.000000000000000111e-01 7.707321840328036373e-05 +2.111111111111111105e-01 8.316552026728771436e-05 +2.222222222222222376e-01 8.636129978080214677e-05 +2.333333333333333370e-01 8.718360874086027276e-05 +2.444444444444444642e-01 8.623107266273786846e-05 +2.555555555555555913e-01 8.387275860560686593e-05 +2.666666666666666630e-01 8.026984295557315521e-05 +2.777777777777777901e-01 7.555489242274063418e-05 +2.888888888888889173e-01 6.991174273120610285e-05 +2.999999999999999889e-01 6.347783338191834092e-05 +3.111111111111111160e-01 5.626067660094025965e-05 +3.222222222222222432e-01 4.828992549022265770e-05 +3.333333333333333703e-01 3.989393478370251709e-05 +3.444444444444444420e-01 3.168527124293345421e-05 +3.555555555555555691e-01 2.408491475986813282e-05 +3.666666666666666963e-01 1.685928263790101785e-05 +3.777777777777777679e-01 9.332803157868673947e-06 +3.888888888888888951e-01 1.221459386476058274e-06 +4.000000000000000222e-01 -6.810527145834715347e-06 +4.111111111111111494e-01 -1.358180240306974179e-05 +4.222222222222222210e-01 -1.838636713962990893e-05 +4.333333333333333481e-01 -2.143290971749338717e-05 +4.444444444444444753e-01 -2.329872060883883496e-05 +4.555555555555555469e-01 -2.416692999531072332e-05 +4.666666666666666741e-01 -2.379516435394250564e-05 +4.777777777777778012e-01 -2.214276017953382861e-05 +4.888888888888889284e-01 -1.974726869296072272e-05 +5.000000000000000000e-01 -1.734254538386544020e-05 +5.111111111111111827e-01 -1.522231256756220318e-05 +5.222222222222222543e-01 -1.313099008534702736e-05 +5.333333333333333259e-01 -1.074976851507047505e-05 +5.444444444444445086e-01 -8.151778532917409610e-06 +5.555555555555555802e-01 -5.748893785602686741e-06 +5.666666666666666519e-01 -3.904220464966048215e-06 +5.777777777777778345e-01 -2.686285556511387382e-06 +5.888888888888889062e-01 -1.952846777477422843e-06 +5.999999999999999778e-01 -1.557977158573065957e-06 +6.111111111111111605e-01 -1.428106984690457702e-06 +6.222222222222222321e-01 -1.491732694493028930e-06 +6.333333333333333037e-01 -1.633848893916421188e-06 +6.444444444444444864e-01 -1.785925925062549406e-06 +6.555555555555555580e-01 -2.051507851562523016e-06 +6.666666666666667407e-01 -2.663863828970359351e-06 +6.777777777777778123e-01 -3.739547050960389822e-06 +6.888888888888888840e-01 -5.076244289743019663e-06 +7.000000000000000666e-01 -6.269804766542887362e-06 +7.111111111111111382e-01 -7.071058511233771672e-06 +7.222222222222222099e-01 -7.575742239734086784e-06 +7.333333333333333925e-01 -8.014752219741598846e-06 +7.444444444444444642e-01 -8.419865022937885718e-06 +7.555555555555555358e-01 -8.600209680439256687e-06 +7.666666666666667185e-01 -8.429971494648805962e-06 +7.777777777777777901e-01 -8.022938266533514274e-06 +7.888888888888889728e-01 -7.566095801927185752e-06 +8.000000000000000444e-01 -7.098279433111552240e-06 +8.111111111111111160e-01 -6.568990729712377425e-06 +8.222222222222222987e-01 -6.037601030251277555e-06 +8.333333333333333703e-01 -5.643254574514988833e-06 +8.444444444444444420e-01 -5.375199879333458125e-06 +8.555555555555556246e-01 -5.056608503029805432e-06 +8.666666666666666963e-01 -4.652311985696004512e-06 +8.777777777777777679e-01 -4.447961326576372164e-06 +8.888888888888889506e-01 -4.748228476792333338e-06 +9.000000000000000222e-01 -5.448030853072816225e-06 +9.111111111111110938e-01 -6.087928154143909244e-06 +9.222222222222222765e-01 -6.355140246645855874e-06 +9.333333333333333481e-01 -6.353499199540407426e-06 +9.444444444444445308e-01 -6.299503166887456157e-06 +9.555555555555556024e-01 -6.133293284942808649e-06 +9.666666666666666741e-01 -5.628094249643691837e-06 +9.777777777777778567e-01 -4.804947720027944289e-06 +9.888888888888889284e-01 -3.992908349549955654e-06 +1.000000000000000000e+00 -3.425824900698486656e-06 +1.011111111111111072e+00 -2.988425260950394039e-06 +1.022222222222222365e+00 -2.502029112748187246e-06 +1.033333333333333437e+00 -2.097731471409650300e-06 +1.044444444444444509e+00 -2.053044079734801498e-06 +1.055555555555555580e+00 -2.274448730253089779e-06 +1.066666666666666652e+00 -2.220509840415812151e-06 +1.077777777777777724e+00 -1.531464435157532165e-06 +1.088888888888889017e+00 -6.373840759172933434e-07 +1.100000000000000089e+00 -5.239489231024034121e-07 diff --git a/test/jointbif/plot.py b/test/jointbif/plot.py new file mode 100644 index 0000000..80b1959 --- /dev/null +++ b/test/jointbif/plot.py @@ -0,0 +1,154 @@ +import numpy as np +import matplotlib.pyplot as plt +import seaborn as sns + + +P = np.loadtxt("bifurcation_results/P_P.out") +Q = np.loadtxt("bifurcation_results/P_Q.out") +A = np.loadtxt("bifurcation_results/P_A.out") +R = np.sqrt(A[:,3]/np.pi) +R0 = 0.83183427183602578*1e-2 +0.000125 +dR = R - R0 + +Rb = np.sqrt(A[:,-1]/np.pi) +R0b = 0.83183427183602578*1e-2 +0.000125 +dRb = Rb - R0b + +P1 = np.loadtxt("bifurcation_results/d1_P.out") +Q1 = np.loadtxt("bifurcation_results/d1_Q.out") +A1 = np.loadtxt("bifurcation_results/d1_A.out") +R1 = np.sqrt(A1[:,3]/np.pi) +R01 = 0.58844816698429576*1e-2 +0.00003 +dR1 = R1 - R01 + + +sns.set_style("white") +sns.set_context("notebook", font_scale=1, rc={"lines.linewidth": 2}) +# plt.locator_params(nbins=5) +fig = plt.figure(1) +fig.clf() + +ax1 = fig.add_subplot(331) +ax1.plot(P[:,0], P[:,3]/1000, 'k') +ax1.set_xlim(8.8, 9.9) +ax1.set_ylim(8, 18) +ax1.set_title("(a) Pressure") + +pad = 5 +ax1.annotate("aorta", xy=(0, 0.5), xytext=(-ax1.yaxis.labelpad - pad, 0), + xycoords=ax1.yaxis.label, textcoords='offset points', + size='large', ha='right', va='center') + +plt.setp( ax1.get_xticklabels(), visible=False) +ax1.set_ylabel("$P$ $($kPa$)$") +ax1.legend(loc="best") + +ax2 = fig.add_subplot(332) +ax2.plot(Q[:,0], Q[:,3]*1e6, 'k') +ax2.set_xlim(8.8, 9.9) +ax2.set_ylim(-24, 80) +ax2.set_title("(b) Flow") +plt.setp( ax2.get_xticklabels(), visible=False) +ax2.set_ylabel("$Q$ $($ml$\cdot$s$^{-1})$") + +ax3 = fig.add_subplot(333) +ax3.plot(A[:,0], dR*1e3, 'k') +ax3.set_xlim(8.8, 9.9) +ax3.set_ylim(-0.05, 0.9) +ax3.set_title("(c) Radius Change") +plt.setp( ax3.get_xticklabels(), visible=False) +ax3.set_ylabel("$\Delta r$ $($mm$)$") + +ax4 = fig.add_subplot(334) +ax4.plot(P[:,0], P[:,-1]/1000, 'k') +ax4.set_xlim(8.8, 9.9) +ax4.set_ylim(7, 19) + +ax4.annotate("bifurcation", xy=(0, 0.5), xytext=(-ax4.yaxis.labelpad - pad, 0), + xycoords=ax4.yaxis.label, textcoords='offset points', + size='large', ha='right', va='center') + +plt.setp( ax4.get_xticklabels(), visible=False) +ax4.set_ylabel("$P$ $($kPa$)$") + +ax5 = fig.add_subplot(335) +ax5.plot(Q[:,0], Q[:,-1]*1e6, 'k') +ax5.set_xlim(8.8, 9.9) +ax5.set_ylim(-20, 70) +plt.setp( ax5.get_xticklabels(), visible=False) +ax5.set_ylabel("$Q$ $($ml$\cdot$s$^{-1})$") + +ax6 = fig.add_subplot(336) +ax6.plot(A[:,0], dRb*1e3, 'k') +ax6.set_xlim(8.8, 9.9) +ax6.set_ylim(-0.1, 0.9) +plt.setp( ax6.get_xticklabels(), visible=False) +ax6.set_ylabel("$\Delta r$ $($mm$)$") + +ax7 = fig.add_subplot(337) +ax7.plot(P1[:,0], P1[:,3]/1000, 'k') +ax7.set_xlim(8.8, 9.9) +ax7.set_ylim(7, 18) + +ax7.annotate("iliac", xy=(0, 0.5), xytext=(-ax7.yaxis.labelpad - pad, 0), + xycoords=ax7.yaxis.label, textcoords='offset points', + size='large', ha='right', va='center') + +ax7.set_ylabel("$P$ $($kPa$)$") +ax7.set_xlabel("$t$ $($s$)$") + +ax8 = fig.add_subplot(338) +ax8.plot(Q1[:,0], Q1[:,3]*1e6, 'k') +ax8.set_xlim(8.8, 9.9) +ax8.set_ylim(-10, 30) +ax8.set_ylabel("$Q$ $($ml$\cdot$s$^{-1})$") +ax8.set_xlabel("$t$ $($s$)$") + +ax9 = fig.add_subplot(339) +ax9.plot(A1[:,0], dR1*1e3, 'k') +ax9.set_xlim(8.8, 9.9) +ax9.set_ylim(-0.07, 0.46) +ax9.set_ylabel("$\Delta r$ $($mm$)$") +ax9.set_xlabel("$t$ $($s$)$") + +sns.despine() + +''' +def quin(t): + return 1e-6 * np.exp(-1e4*(t-0.05)**2) + +t = np.linspace(0., 2., 1000) +quint = quin(t) + +quint += 1e-10 + +data = np.vstack((t, quint)) + +np.savetxt("tutorial_inlet.dat", np.transpose(data)) + +plt.figure(2) +plt.clf() +plt.plot(t, quint) + + +PI=3.141592653589793 +period T=1.1 +time t +inlet flow rate in ml/s + +T = 1.1 +t = np.linspace(0., T, 100) +PI = np.pi +inflow=10e5*(7.9853e-06+2.6617e-05*np.sin(2*PI*t/T+0.29498)+2.3616e-05*np.sin(4*PI*t/T-1.1403)-1.9016e-05*np.sin(6*PI*t/T+0.40435)-8.5899e-06*np.sin(8*PI*t/T-1.1892)-2.436e-06*np.sin(10*PI*t/T-1.4918)+1.4905e-06*np.sin(12*PI*t/T+1.0536)+1.3581e-06*np.sin(14*PI*t/T-0.47666)-6.3031e-07*np.sin(16*PI*t/T+0.93768)-4.5335e-07*np.sin(18*PI*t/T-0.79472)-4.5184e-07*np.sin(20*PI*t/T-1.4095)-5.6583e-07*np.sin(22*PI*t/T-1.3629)+4.9522e-07*np.sin(24*PI*t/T+0.52495)+1.3049e-07*np.sin(26*PI*t/T-0.97261)-4.1072e-08*np.sin(28*PI*t/T-0.15685)-2.4182e-07*np.sin(30*PI*t/T-1.4052)-6.6217e-08*np.sin(32*PI*t/T-1.3785)-1.5511e-07*np.sin(34*PI*t/T-1.2927)+2.2149e-07*np.sin(36*PI*t/T+0.68178)+6.7621e-08*np.sin(38*PI*t/T-0.98825)+1.0973e-07*np.sin(40*PI*t/T+1.4327)-2.5559e-08*np.sin(42*PI*t/T-1.2372)-3.5079e-08*np.sin(44*PI*t/T+0.2328)) +data = np.vstack((t, inflow*1e-6)) + +np.savetxt("tutorial_inlet.dat", np.transpose(data)) + + +plt.figure(2) +plt.clf() +plt.plot(t, inflow) +''' + +plt.draw() +plt.show() diff --git a/test/jointcat/conjunction.yml b/test/jointcat/conjunction.yml index 30abaa6..7700a88 100644 --- a/test/jointcat/conjunction.yml +++ b/test/jointcat/conjunction.yml @@ -16,10 +16,10 @@ network: sn: 1 # proximal (source) node tn: 2 # distal (target) node - L: 12.14e-2 # length (m) + L: 4.14e-2 # length (m) R0: 9.87e-3 # proximal lumen radius (m) - #Rc: 1.48e-3 # catheter radius + Rc: 7.48e-3 # catheter radius E: 700.0e3 # Young's modulus (Pa) h0: 0.82e-3 @@ -33,7 +33,7 @@ network: sn: 2 tn: 3 - L: 12.14e-2 # length (m) + L: 4.14e-2 # length (m) R0: 9.87e-3 # proximal lumen radius (m) diff --git a/test/jointcat/script.jl b/test/jointcat/script.jl index a6cba03..f4a52d3 100644 --- a/test/jointcat/script.jl +++ b/test/jointcat/script.jl @@ -1,7 +1,7 @@ using Revise using openBF -cd("/home/ivan/repos/openBF/test/single-artery/") -data = openBF.loadSimulationFiles("single-artery.yml") +cd("/home/ivan/Dropbox/Postdoc/Cardiovascular/openBF/oBF_admin/openBF/test/jointcat/") +data = openBF.loadSimulationFiles("conjunction.yml") blood = openBF.buildBlood(data["blood"]) jump = data["solver"]["jump"] vessels, edges = openBF.buildArterialNetwork(data["network"], blood, jump) From 7d281befe75b71c0fb0674ed8c2425eae1ec0ee3 Mon Sep 17 00:00:00 2001 From: ivan Date: Thu, 16 Jan 2020 00:13:01 +0000 Subject: [PATCH 4/4] modified c viscous term --- src/boundary_conditions.jl | 2 +- src/initialise.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/boundary_conditions.jl b/src/boundary_conditions.jl index 683ba94..72eea08 100644 --- a/src/boundary_conditions.jl +++ b/src/boundary_conditions.jl @@ -137,7 +137,7 @@ function setOutletBC(dt :: Float64, v :: Vessel) outletCompatibility(dt, v) elseif v.outlet == "wk3" threeElementWindkessel(dt, v) - end + end end diff --git a/src/initialise.jl b/src/initialise.jl index 5b61f5d..af762e3 100644 --- a/src/initialise.jl +++ b/src/initialise.jl @@ -586,12 +586,12 @@ assumed equal to `9` (plug-flow). function computeViscousTerm(vessel_data :: Dict{Any,Any}, blood :: Blood) if haskey(vessel_data,"Rc") # MODIFIED THIS LINE - return 2*blood.mu*blood.rho_inv*sqrt(15) + return 2*blood.mu*blood.rho_inv*sqrt(30) elseif haskey(vessel_data, "gamma_profile") gamma_profile = vessel_data["gamma_profile"] else gamma_profile = 9 - return 2*(gamma_profile + 2)*pi*blood.mu*blood.rho_inv + return 2*(gamma_profile + 2)*blood.mu*blood.rho_inv end end