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

Initial solution for mesons

This commit is contained in:
Nuwan Yapa 2025-01-30 18:36:29 -05:00
parent 174dc40f20
commit c9ec9c6203
2 changed files with 18 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
mesons.jl
View File

@ -44,10 +44,11 @@ end
ρ_sp, ρ_vp are the scalar and vector proton densities as arrays, ρ_sp, ρ_vp are the scalar and vector proton densities as arrays,
ρ_sn, ρ_vn are the scalar and vector neutron densities as arrays, ρ_sn, ρ_vn are the scalar and vector neutron densities as arrays,
Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as as arrays, Φ0, W0, B0, A0 are the mean-field potentials (couplings included) in MeV as as arrays,
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!).
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) function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, iterations=3; initial_sol=(zeros(1 + divs) for _ in 1:4))
(Φ0, W0, B0, A0) = (zeros(1 + divs) for _ in 1:4) (Φ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)
# A0 doesn't need iterations # A0 doesn't need iterations

26
test/Pb208_meson_flds.jl
View File

@ -1,6 +1,17 @@
using DelimitedFiles, Interpolations, Plots using DelimitedFiles, Interpolations, Plots
include("../mesons.jl") include("../mesons.jl")
# test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x)
bench_data = readdlm("test/Pb208FldsFSUGarnet.csv")
xs_bench = bench_data[:, 1]
Φ0_bench = bench_data[:, 2]
W0_bench = bench_data[:, 3]
B0_bench = bench_data[:, 4]
A0_bench = bench_data[:, 5]
plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["Φ0 benchmark" "W0 benchmark" "B0 benchmark" "A0 benchmark"])
# test data generated from Hartree.f # test data generated from Hartree.f
# format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p) # format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p)
test_data = readdlm("test/Pb208DensFSUGarnet.csv") test_data = readdlm("test/Pb208DensFSUGarnet.csv")
@ -13,18 +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) (Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, r_max, divs, 1; initial_sol=(Φ0_bench, W0_bench, B0_bench, A0_bench))
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)")
# test data generated from Hartree.f
# format: x S(x) V(x) R(x) A(x)
test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
xs_bench = test_data[:, 1]
Φ0_bench = test_data[:, 2]
W0_bench = test_data[:, 3]
B0_bench = test_data[:, 4]
A0_bench = test_data[:, 5]
plot!(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["Φ0 benchmark" "W0 benchmark" "B0 benchmark" "A0 benchmark"])