yapanuwan
/
NuclearRMF
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

KG instability fixed

This commit is contained in:
Nuwan Yapa 2025-02-05 13:49:40 -05:00
parent 32dc165f6c
commit 0c5b906e85
2 changed files with 4 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
mesons.jl
View File

@ -47,7 +47,7 @@ end
r is the radius in fm, r is the radius in fm,
An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!). An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!).
Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)" Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3; initial_sol=(zeros(1 + divs) for _ in 1:4)) function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4))
(Φ0, W0, B0, A0) = initial_sol (Φ0, W0, B0, A0) = initial_sol
(src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4) (src_Φ0, src_W0, src_B0, src_A0) = (zeros(1 + divs) for _ in 1:4)
@ -60,8 +60,8 @@ function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3; in
@. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (2B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn) @. src_W0 = g2_v * ((ζ/6) * (W0/ħc)^3 + 2Λv * (2B0/ħc)^2 * (W0/ħc)) - g2_v * (ρ_vp + ρ_vn)
@. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2 @. src_B0 = (2Λv * g2_ρ * (W0/ħc)^2 * (2B0/ħc) - g2_ρ/2 * (ρ_vp - ρ_vn)) / 2
Φ0 .= solveKG(m_s, src_Φ0, r_max) Φ0 .= solveKG(m_s, src_Φ0, r_max)
W0 .= solveKG(m_ω, src_W0, r_max) W0 .= (solveKG(m_ω, src_W0, r_max) .+ W0) ./ 2 # to supress oscillation
B0 .= solveKG(m_ρ, src_B0, r_max) B0 .= (solveKG(m_ρ, src_B0, r_max) .+ B0) ./ 2 # to supress oscillation
end end
return (Φ0, W0, B0, A0) return (Φ0, W0, B0, A0)

2
test/Pb208_meson_flds.jl
View File

@ -24,7 +24,7 @@ xs = test_data[:, 1]
r_max = maximum(xs) r_max = maximum(xs)
divs = length(xs) - 1 divs = length(xs) - 1
(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 1; initial_sol=(Φ0_bench, W0_bench, B0_bench, A0_bench)) (Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 200) # many iterations needed without an initial solution
plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"]) plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"])
xlabel!("r (fm)") xlabel!("r (fm)")