diff --git a/mesons.jl b/mesons.jl index 243aa56..79e1903 100644 --- a/mesons.jl +++ b/mesons.jl @@ -38,14 +38,13 @@ end "Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where divs is the number of mesh divisions so the arrays are of length (1+divs), - ρ_sp, ρ_vp are the scalar and vector proton 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, r is the radius in fm, - An initial guess initial_sol=(Φ0, W0, B0, A0) can be provided to speed up convergence (permuting!). + the inital solutions are read from s and the final solutions are saved in-place. Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)" -function solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s::system, iterations=10; initial_sol=(zeros(1 + divs) for _ in 1:4)) - (Φ0, W0, B0, A0) = initial_sol +function solveMesonFields!(s::system, iterations=10) + (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0) + (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn) + (src_Φ0, src_W0, src_B0, src_A0) = (zero_array(s) for _ in 1:4) # A0 doesn't need iterations diff --git a/nucleons.jl b/nucleons.jl index dc4fc2e..a1298cc 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -27,7 +27,9 @@ end divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix, shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max, the other parameters are the same from dirac!(...)." -function solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, s::system; shooting=true, normalize=true) +function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true) + (Φ0, W0, B0, A0) = fields_interp(s) + if shooting @assert s.divs % 2 == 0 "divs must be an even number when shooting=true" prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2)) @@ -62,8 +64,9 @@ end "Returns a function that solves the Dirac equation and returns g(r=r_max) where r_max is the outer boundary in fm, the other parameters are the same from dirac!(...)." -function boundaryValueFunc(κ, p, Φ0, W0, B0, A0, s::system; dtype=Float64, algo=Tsit5()) +function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5()) prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max)) + (Φ0, W0, B0, A0) = fields_interp(s) func(E) = solve(prob, algo, p=(κ, p, E, Φ0, W0, B0, A0), saveat=[s.r_max], save_idxs=[1])[1, 1] return func end @@ -71,21 +74,21 @@ end "Find all bound energies between E_min (=0) and E_max (=mass) where tol_digits determines the precision for root finding and the threshold for identifying duplicates, the other parameters are the same from dirac!(...)." -function findEs(κ, p, Φ0, W0, B0, A0, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5) - func = boundaryValueFunc(κ, p, Φ0, W0, B0, A0, s) +function findEs(κ, p, s::system, E_min=0, E_max=(p ? M_p : M_n), tol_digits=5) + func = boundaryValueFunc(κ, p, s) Es = find_zeros(func, (E_min, E_max); xatol=1/10^tol_digits) return unique(E -> round(E; digits=tol_digits), Es) end "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where the other parameters are defined above" -function findAllOrbitals(p, Φ0, W0, B0, A0, s::system, E_min=0, E_max=(p ? M_p : M_n)) +function findAllOrbitals(p, s::system, E_min=0, E_max=(p ? M_p : M_n)) κs = Int[] Es = Float64[] # start from κ=-1 and go both up and down for direction in [-1, 1] for κ in direction * (1:100) # cutoff is 100 - new_Es = findEs(κ, p, Φ0, W0, B0, A0, s, E_min, E_max) + new_Es = findEs(κ, p, s, E_min, E_max) if isempty(new_Es); break; end append!(Es, new_Es) append!(κs, fill(κ, length(new_Es))) @@ -120,11 +123,11 @@ l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0 "Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ, the other parameters are defined above" -function calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, s::system) +function calculateNucleonDensity(κs, Es, occs, p, s::system) ρr2 = zeros(2, s.divs + 1) # ρ×r² values for (κ, E, occ) in zip(κs, Es, occs) - wf = solveNucleonWf(κ, p, E, Φ0, W0, B0, A0, s; shooting=true, normalize=true) + wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true) wf2 = wf .* wf ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number end @@ -141,8 +144,8 @@ end "Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where the other parameters are defined above" -function solveNucleonDensity(p, Φ0, W0, B0, A0, s::system, E_min=800, E_max=939) - κs, Es = findAllOrbitals(p, Φ0, W0, B0, A0, s, E_min, E_max) +function solveNucleonDensity(p, s::system, E_min=800, E_max=939) + κs, Es = findAllOrbitals(p, s, E_min, E_max) occs = fillNucleons(Z_or_N(s, p), κs, Es) - return calculateNucleonDensity(κs, Es, occs, p, Φ0, W0, B0, A0, s) + return calculateNucleonDensity(κs, Es, occs, p, s) end diff --git a/system.jl b/system.jl index 4490727..c3a2531 100644 --- a/system.jl +++ b/system.jl @@ -18,11 +18,11 @@ mutable struct system ρ_sn::Vector{Float64} ρ_vn::Vector{Float64} - "Dummy struct to define the mesh" - system(r_max, divs) = new(0, 0, r_max, divs, [Float64[] for _ in 1:8]...) - "Initialize an unsolved system" system(Z, N, r_max, divs) = new(Z, N, r_max, divs, [zeros(1 + divs) for _ in 1:8]...) + + "Dummy struct to define the mesh" + system(r_max, divs) = system(0, 0, r_max, divs) end "Get mass number of nucleus" @@ -73,18 +73,17 @@ function solve_system(s::system; reinitialize_densities=true, monitor_print=true while true # mesons - @time "Meson fields" solveMesonWfs(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s, isnan(E_total_previous) ? 50 : 5; initial_sol = (s.Φ0, s.W0, s.B0, s.A0)) - (S_interp, V_interp, R_interp, A_interp) = fields_interp(s) + @time "Meson fields" solveMesonFields!(s, isnan(E_total_previous) ? 50 : 5) # protons - @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, S_interp, V_interp, R_interp, A_interp, s) + @time "Proton spectrum" (κs_p, Es_p) = findAllOrbitals(true, s) occs_p = fillNucleons(s.Z, κs_p, Es_p) - @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, S_interp, V_interp, R_interp, A_interp, s) + @time "Proton densities" (s.ρ_sp, s.ρ_vp) = calculateNucleonDensity(κs_p, Es_p, occs_p, true, s) # neutrons - @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, S_interp, V_interp, R_interp, A_interp, s) + @time "Neutron spectrum" (κs_n, Es_n) = findAllOrbitals(false, s) occs_n = fillNucleons(s.N, κs_n, Es_n) - @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, S_interp, V_interp, R_interp, A_interp, s) + @time "Neutron densities" (s.ρ_sn, s.ρ_vn) = calculateNucleonDensity(κs_n, Es_n, occs_n, false, s) if monitor_plot for s in p.series_list diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl index ffbf48f..c438106 100644 --- a/test/Pb208_meson_flds.jl +++ b/test/Pb208_meson_flds.jl @@ -16,16 +16,17 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=[" # format: x Rhos(n) Rhov(n) Rhot(n) Rhos(p) Rhov(p) Rhot(p) test_data = readdlm("test/Pb208DensFSUGarnet.csv") xs = test_data[:, 1] -ρ_sn = test_data[:, 2] -ρ_vn = test_data[:, 3] -ρ_sp = test_data[:, 5] -ρ_vp = test_data[:, 6] r_max = maximum(xs) divs = length(xs) - 1 s = system(r_max, divs) -(Φ0, W0, B0, A0) = solveMesonWfs(ρ_sp, ρ_vp, ρ_sn, ρ_vn, s, 200) # many iterations needed without an initial solution +s.ρ_sn = test_data[:, 2] +s.ρ_vn = test_data[:, 3] +s.ρ_sp = test_data[:, 5] +s.ρ_vp = test_data[:, 6] + +(Φ0, W0, B0, A0) = solveMesonFields!(s, 200) # many iterations needed without an initial solution plot!(xs, hcat(Φ0, W0, B0, A0), layout=4, label=["Φ0" "W0" "B0" "A0"]) xlabel!("r (fm)") \ No newline at end of file diff --git a/test/Pb208_nucleon_basic.jl b/test/Pb208_nucleon_basic.jl index 6349a39..a5f2edf 100644 --- a/test/Pb208_nucleon_basic.jl +++ b/test/Pb208_nucleon_basic.jl @@ -21,17 +21,15 @@ s.W0 = Vs s.B0 = Rs s.A0 = As -(S_interp, V_interp, R_interp, A_interp) = fields_interp(s) - E_min = 850 E_max = 939 -boundEs = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max) +boundEs = findEs(κ, p, s, E_min, E_max) binding_Es = (p ? M_p : M_n) .- boundEs println("binding energies = $binding_Es") -func = boundaryValueFunc(κ, p, S_interp, V_interp, R_interp, A_interp, s) +func = boundaryValueFunc(κ, p, s) Es = collect(E_min:0.5:E_max) boundaryVals = [func(E)^2 for E in Es] diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl index 6ffb7a9..9c1c277 100644 --- a/test/Pb208_nucleon_dens.jl +++ b/test/Pb208_nucleon_dens.jl @@ -21,14 +21,12 @@ s.W0 = Vs s.B0 = Rs s.A0 = As -(S_interp, V_interp, R_interp, A_interp) = fields_interp(s) - -(ρ_sp, ρ_vp) = solveNucleonDensity(true, S_interp, V_interp, R_interp, A_interp, s) +(ρ_sp, ρ_vp) = solveNucleonDensity(true, s) p_sp = plot(rs(s), ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated") p_vp = plot(rs(s), ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated") -(ρ_sn, ρ_vn) = solveNucleonDensity(false, S_interp, V_interp, R_interp, A_interp, s) +(ρ_sn, ρ_vn) = solveNucleonDensity(false, s) p_sn = plot(rs(s), ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated") p_vn = plot(rs(s), ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated") diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl index ed20529..c12ea6d 100644 --- a/test/Pb208_nucleon_spectrum.jl +++ b/test/Pb208_nucleon_spectrum.jl @@ -22,12 +22,10 @@ s.W0 = Vs s.B0 = Rs s.A0 = As -(S_interp, V_interp, R_interp, A_interp) = fields_interp(s) - E_min = 800 E_max = 939 -κs, Es = findAllOrbitals(p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max) +κs, Es = findAllOrbitals(p, s, E_min, E_max) Ebinds = (p ? M_p : M_n) .- Es occ = fillNucleons(Z_or_N(s, p), κs, Es) diff --git a/test/Pb208_nucleon_wf.jl b/test/Pb208_nucleon_wf.jl index 59543f4..7898949 100644 --- a/test/Pb208_nucleon_wf.jl +++ b/test/Pb208_nucleon_wf.jl @@ -21,15 +21,13 @@ s.W0 = Vs s.B0 = Rs s.A0 = As -(S_interp, V_interp, R_interp, A_interp) = fields_interp(s) - E_min = 800 E_max = 939 -groundE = findEs(κ, p, S_interp, V_interp, R_interp, A_interp, s, E_min, E_max) |> minimum +groundE = findEs(κ, p, s, E_min, E_max) |> minimum println("ground state E = $groundE") -wf = solveNucleonWf(κ, p, groundE, S_interp, V_interp, R_interp, A_interp, s) +wf = solveNucleonWf(κ, p, groundE, s) gs = wf[1, :] fs = wf[2, :]