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

Now supports predefined mesh

This commit is contained in:
Nuwan Yapa 2025-03-12 19:02:50 -04:00
parent 33c8b3f81d
commit 2dc59a00d2
6 changed files with 63 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
NuclearRMF.jl
View File

@ -11,7 +11,7 @@ Woods_Saxon(r::Float64; R::Float64=7.0, a::Float64=0.5) = -1 / (8π * a^3 * reli
"Run the full program by self-consistent solution of nucleon and meson densities" "Run the full program by self-consistent solution of nucleon and meson densities"
function solve_system!(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false) function solve_system!(s::system; reinitialize_densities=true, monitor_print=true, monitor_plot=false)
if reinitialize_densities if reinitialize_densities
dens_guess = Woods_Saxon.(rs(s)) dens_guess = Woods_Saxon.(s.r_mesh.r)
@. s.ρ_sp = s.Z * dens_guess @. s.ρ_sp = s.Z * dens_guess
@. s.ρ_vp = s.Z * dens_guess @. s.ρ_vp = s.Z * dens_guess
@. s.ρ_sn = s.N * dens_guess @. s.ρ_sn = s.N * dens_guess
@ -33,7 +33,7 @@ function solve_system!(s::system; reinitialize_densities=true, monitor_print=tru
for s in p.series_list for s in p.series_list
s.plotattributes[:linecolor] = :gray s.plotattributes[:linecolor] = :gray
end end
plot!(p, rs(s), hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red) plot!(p, s.r_mesh.r, hcat(s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn, s.Φ0, s.W0, s.B0, s.A0), linecolor=:red)
display(p) display(p)
end end

11
mesons.jl
View File

@ -25,16 +25,15 @@ greensFunctionP(r, rp) = -1 / (max(r, rp) + r_reg)
"Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where "Solve the Klein-Gordon equation (or Poisson's equation when m=0) and return an array in MeV for a source array given in fm⁻³ where
m is the mass of the meson in MeV/c2." m is the mass of the meson in MeV/c2."
function solveKG(m, source, s::system) function solveKG(m, source, s::system)
int_measure = ħc .* Δr(s) .* rs(s) .^ 2 int_measure = ħc .* s.r_mesh.w .* (s.r_mesh.r .^ 2)
greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp) greensFunction = m == 0 ? greensFunctionP : (r, rp) -> greensFunctionKG(m / ħc, r, rp)
greenMat = greensFunction.(rs(s), transpose(rs(s))) greenMat = greensFunction.(s.r_mesh.r, transpose(s.r_mesh.r))
return greenMat * (int_measure .* source) return greenMat * (int_measure .* source)
end end
"Iteratively solve meson equations and return the wave functions u(r)=[Φ0(r), W0(r), B0(r), A0(r)] where "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),
r is the radius in fm, r is the radius in fm,
the inital solutions are read from s and the final solutions are saved in-place. 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)" Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
@ -69,13 +68,13 @@ end
"Calculate the total energy associated with meson fields" "Calculate the total energy associated with meson fields"
function meson_E(s::system) function meson_E(s::system)
int = 0.0 int = 0.0
for (r, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(rs(s), s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn) for (r, w, Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn) in zip(s.r_mesh.r, s.r_mesh.w, s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4) E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4 E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn) E_ρ = -(1/4) * (2B0/ħc) * (ρ_vp - ρ_vn)
E_γ = -(1/2) * (A0/ħc) * ρ_vp E_γ = -(1/2) * (A0/ħc) * ρ_vp
E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2 E_ωρ = Λv * (W0/ħc)^2 * (2B0/ħc)^2
int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2 int += (E_σ + E_ω + E_ρ + E_γ + E_ωρ) * r^2 * w
end end
return 4π * int * Δr(s) * ħc return 4π * int * ħc
end end

63
nucleons.jl
View File

@ -33,44 +33,45 @@ function optimized_dirac_potentials(p::Bool, s::system)
@. f1s = M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0 @. f1s = M - s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
@. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0 @. f2s = -M + s.Φ0 - s.W0 - (p - 0.5) * 2 * s.B0 - p * s.A0
f1 = linear_interpolation(rs(s), f1s) f1 = linear_interpolation(s.r_mesh.r, f1s)
f2 = linear_interpolation(rs(s), f2s) f2 = linear_interpolation(s.r_mesh.r, f2s)
return (f1, f2) return (f1, f2)
end end
"Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where "Solve the Dirac equation and return the wave function u(r)=[g(r), f(r)] where
divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix, the solution would be returned as a 2×mesh_size matrix,
shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max, 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!(...)." the other parameters are the same from dirac!(...)."
function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Tsit5()) function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Tsit5())
(f1, f2) = optimized_dirac_potentials(p, s) (f1, f2) = optimized_dirac_potentials(p, s)
if shooting # partitioning
@assert s.divs % 2 == 0 "divs must be an even number when shooting=true" mid_idx = length(s.r_mesh.r) ÷ 2
prob = ODEProblem(dirac!, [0, 1], (s.r_max, s.r_max / 2)) r_mid = s.r_mesh.r[mid_idx]
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) left_r = s.r_mesh.r[1:mid_idx]
wf_right = reverse(hcat(sol.u...); dims=2) right_r = s.r_mesh.r[mid_idx:end]
next_r_max = s.r_max / 2 # for the next step
else
next_r_max = s.r_max
end
prob = ODEProblem(dirac!, [0, 1], (0, next_r_max)) # left partition
sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) prob = ODEProblem(dirac!, [0, 1], (0, r_mid))
wf = hcat(sol.u...) sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r)
wf_left = hcat(sol.u...)
if shooting # join two segments # right partition
rescale_factor_g = wf[1, end] / wf_right[1, 1] prob = ODEProblem(dirac!, [0, 1], (s.r_mesh.r_max, r_mid))
rescale_factor_f = wf[2, end] / wf_right[2, 1] sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r)
@assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition" wf_right = reverse(hcat(sol.u...); dims=2)
isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2) # join two segments
wf = hcat(wf[:, 1:(end - 1)], wf_right_rescaled) rescale_factor_g = wf_left[1, end] / wf_right[1, 1]
end rescale_factor_f = wf_left[2, end] / wf_right[2, 1]
@assert isfinite(rescale_factor_g) && isfinite(rescale_factor_f) "Cannot rescale the right partition"
isapprox(rescale_factor_g, rescale_factor_f; rtol=0.03) || @warn "Discontinuity between two partitions"
wf_right_rescaled = wf_right .* ((rescale_factor_g + rescale_factor_f) / 2)
wf = hcat(wf_left[:, 1:(end - 1)], wf_right_rescaled)
if normalize if normalize
norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum norm = sum(wf .* wf .* transpose(s.r_mesh.w)) # integration by Reimann sum
wf = wf ./ sqrt(norm) wf = wf ./ sqrt(norm)
end end
@ -80,10 +81,10 @@ end
"Returns a function that solves the Dirac equation in two partitions and returns "Returns a function that solves the Dirac equation in two partitions and returns
the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)], the determinant of [g_left(r) g_right(r); f_left(r) f_right(r)],
where is r is in fm." where is r is in fm."
function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Tsit5()) function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_mesh.r_max/2, algo=Tsit5())
(f1, f2) = optimized_dirac_potentials(p, s) (f1, f2) = optimized_dirac_potentials(p, s)
prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r)) prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r))
prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_max, r)) prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_mesh.r_max, r))
function func(E) function func(E)
u_left = solve(prob_left, algo, p=(κ, E, f1, f2), saveat=[r]) u_left = solve(prob_left, algo, p=(κ, E, f1, f2), saveat=[r])
u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r]) u_right = solve(prob_right, algo, p=(κ, E, f1, f2), saveat=[r])
@ -141,22 +142,22 @@ j_κ(κ::Int) = abs(κ) - 1/2
"Orbital angular momentum l for a given κ value" "Orbital angular momentum l for a given κ value"
l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0 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 "Calculate scalar and vector densities of a nucleon species evaluated at mesh points and returns them as vectors (ρ_s, ρ_v) where
the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ, the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
the other parameters are defined above" the other parameters are defined above"
function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}} function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vector{Float64}}
spectrum = p ? s.p_spectrum : s.n_spectrum spectrum = p ? s.p_spectrum : s.n_spectrum
(κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ) (κs, Es, occs) = (spectrum.κ, spectrum.E, spectrum.occ)
ρr2 = zeros(2, s.divs + 1) # ρ×r² values ρr2 = zeros(2, length(s.r_mesh.r)) # ρ×r² values
for (κ, E, occ) in zip(κs, Es, occs) for (κ, E, occ) in zip(κs, Es, occs)
wf = solveNucleonWf(κ, p, E, s; shooting=true, normalize=true) wf = solveNucleonWf(κ, p, E, s; normalize=true)
wf2 = wf .* wf wf2 = wf .* wf
ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number
end end
r2s = rs(s).^2 r2s = s.r_mesh.r .^ 2
ρ = ρr2 ./ transpose(r2s) ρ = ρr2 ./ transpose(r2s)
ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0 ρ[:, 1] .= ρ[:, 2] # dirty fix for NaN at r=0

27
system.jl
View File

@ -10,13 +10,26 @@ end
"Initializes an unfilled spectrum" "Initializes an unfilled spectrum"
unfilled_spectrum() = spectrum(Int[], Float64[], Int[]) unfilled_spectrum() = spectrum(Int[], Float64[], Int[])
"Defines a mesh"
struct mesh
r_max::Float64
r::Vector{Float64}
w::Vector{Float64}
end
"Create a uniform mesh"
function uniform_mesh(r_max::Float64, divs::Int)::mesh
r = range(0, r_max, length=divs+1) |> collect
w = fill(r_max / divs, divs+1)
return mesh(r_max, r, w)
end
"Defines a nuclear system containing relevant parameters and meson/nucleon densities" "Defines a nuclear system containing relevant parameters and meson/nucleon densities"
mutable struct system mutable struct system
Z::Int Z::Int
N::Int N::Int
r_max::Float64 r_mesh::mesh
divs::Int
p_spectrum::spectrum p_spectrum::spectrum
n_spectrum::spectrum n_spectrum::spectrum
@ -32,7 +45,7 @@ mutable struct system
ρ_vn::Vector{Float64} ρ_vn::Vector{Float64}
"Initialize an unsolved system" "Initialize an unsolved system"
system(Z, N, r_max, divs) = new(Z, N, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...) system(Z, N, r_max, divs) = new(Z, N, uniform_mesh(r_max, divs), unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)
"Dummy struct to define the mesh" "Dummy struct to define the mesh"
system(r_max, divs) = system(0, 0, r_max, divs) system(r_max, divs) = system(0, 0, r_max, divs)
@ -41,14 +54,8 @@ end
"Get mass number of nucleus" "Get mass number of nucleus"
A(s::system)::Int = s.Z + s.N A(s::system)::Int = s.Z + s.N
"Get r values in the mesh"
rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)
"Get Δr value for the mesh"
Δr(s::system)::Float64 = s.r_max / s.divs
"Get the number of protons or neutrons in the system" "Get the number of protons or neutrons in the system"
Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N
"Create an empty array for the size of the mesh" "Create an empty array for the size of the mesh"
zero_array(s::system) = zeros(1 + s.divs) zero_array(s::system) = zeros(length(s.r_mesh.r))

8
test/Pb208_nucleon_dens.jl
View File

@ -23,13 +23,13 @@ s.A0 = As
solveNucleonDensity!(true, s) solveNucleonDensity!(true, s)
p_sp = plot(rs(s), s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated") p_sp = plot(s.r_mesh.r, s.ρ_sp, xlabel="r (fm)", label="ρₛₚ(r) calculated")
p_vp = plot(rs(s), s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated") p_vp = plot(s.r_mesh.r, s.ρ_vp, xlabel="r (fm)", label="ρᵥₚ(r) calculated")
solveNucleonDensity!(false, s) solveNucleonDensity!(false, s)
p_sn = plot(rs(s), s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated") p_sn = plot(s.r_mesh.r, s.ρ_sn, xlabel="r (fm)", label="ρₛₙ(r) calculated")
p_vn = plot(rs(s), s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated") p_vn = plot(s.r_mesh.r, s.ρ_vn, xlabel="r (fm)", label="ρᵥₙ(r) calculated")
# benchmark data generated from Hartree.f # benchmark 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)

4
test/Pb208_nucleon_wf.jl
View File

@ -28,6 +28,6 @@ wf = solveNucleonWf(κ, p, groundE, s)
gs = wf[1, :] gs = wf[1, :]
fs = wf[2, :] fs = wf[2, :]
plot(rs(s), gs, label="g(r)") plot(s.r_mesh.r, gs, label="g(r)")
plot!(rs(s), fs, label="f(r)") plot!(s.r_mesh.r, fs, label="f(r)")
xlabel!("r (fm)") xlabel!("r (fm)")