From 9f8cea12a3296236cd05f4c1be9d951b7ef11ea4 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Fri, 4 Apr 2025 18:36:01 -0400 Subject: [PATCH 1/6] Asymptotic boundary condition --- nucleons.jl | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/nucleons.jl b/nucleons.jl index 4b92e20..c70a9be 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -39,6 +39,18 @@ function optimized_dirac_potentials(p::Bool, s::system) return (f1, f2) end +"Approximate boundary condition for u(r)=[g(r), f(r)] at r -> ∞ where + κ is the generalized angular momentum, + p is true for proton and false for neutron, + E is the energy in MeV, + r is the radius in fm." +function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64) + M = p ? M_p : M_n + g = 1.0 + f = ħc / (E + M) * (-√(M^2 - E^2) + κ/r) * g + return [g, f] +end + "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, shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max, @@ -48,7 +60,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true 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)) + prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, s.r_max / 2)) sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) wf_right = reverse(hcat(sol.u...); dims=2) next_r_max = s.r_max / 2 # for the next step @@ -83,8 +95,8 @@ end function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Tsit5()) (f1, f2) = optimized_dirac_potentials(p, s) prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r)) - prob_right = ODEProblem(dirac!, [0.0, 1.0], (s.r_max, r)) function func(E) + prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, 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]) return u_left[1, 1] * u_right[2, 1] - u_right[1, 1] * u_left[2, 1] From a143edb49630f52adfcb7b44bae365970c07579b Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Mon, 7 Apr 2025 19:24:42 -0400 Subject: [PATCH 2/6] Elegant method to attach two partitions --- nucleons.jl | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/nucleons.jl b/nucleons.jl index c70a9be..d465af4 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -1,4 +1,4 @@ -using DifferentialEquations, Interpolations +using LinearAlgebra, DifferentialEquations, Interpolations include("bisection.jl") include("common.jl") include("system.jl") @@ -73,17 +73,20 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true wf = hcat(sol.u...) if shooting # join two segments - rescale_factor_g = wf[1, end] / wf_right[1, 1] - rescale_factor_f = wf[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[:, 1:(end - 1)], wf_right_rescaled) + u1 = wf[:, end] + u2 = wf_right[:, 1] + if norm(u2) < 1e-10 + @warn "Right partition too small to rescale, setting to zero" + wf_right .= 0.0 + else + proj = only(u1' * u2) / norm(u2)^2 + wf_right .*= proj + end + wf = hcat(wf[:, 1:(end - 1)], wf_right) end if normalize - norm = sum(wf .* wf) * Δr(s) # integration by Reimann sum - wf = wf ./ sqrt(norm) + wf ./= norm(wf) * Δr(s) # integration by Reimann sum end return wf From d3b2ec4dd8d7441200e5f4b95b5858d7e63ecb22 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Mon, 7 Apr 2025 19:38:23 -0400 Subject: [PATCH 3/6] Integration bug fixed --- nucleons.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nucleons.jl b/nucleons.jl index d465af4..4796bbf 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -86,7 +86,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true end if normalize - wf ./= norm(wf) * Δr(s) # integration by Reimann sum + wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum end return wf From 268ce12b6cc7415af09f5b1e948832b9a9959748 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Mon, 7 Apr 2025 19:40:31 -0400 Subject: [PATCH 4/6] Tsit5 -> Vern9 for better accuary --- nucleons.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nucleons.jl b/nucleons.jl index 4796bbf..7fbc194 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -55,7 +55,7 @@ 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::Bool, E, s::system; shooting=true, normalize=true, algo=Tsit5()) +function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Vern9()) (f1, f2) = optimized_dirac_potentials(p, s) if shooting @@ -95,7 +95,7 @@ end "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)], 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_max/2, algo=Vern9()) (f1, f2) = optimized_dirac_potentials(p, s) prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r)) function func(E) From 897fdf921496d2ef2835092f71c9436f349f91de Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Tue, 8 Apr 2025 12:29:01 -0400 Subject: [PATCH 5/6] Changed initial boundary condition --- nucleons.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/nucleons.jl b/nucleons.jl index 7fbc194..3c3f81d 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -51,6 +51,9 @@ function asymp_BC(κ::Int, p::Bool, E::Float64, r::Float64) return [g, f] end +"Initial boundary condition for u(r)=[g(r), f(r)] at r=0" +init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]? + "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, shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max, @@ -68,7 +71,7 @@ function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true next_r_max = s.r_max end - prob = ODEProblem(dirac!, [0, 1], (0, next_r_max)) + prob = ODEProblem(dirac!, init_BC(), (0, next_r_max)) sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) wf = hcat(sol.u...) @@ -97,7 +100,7 @@ end where is r is in fm." function determinantFunc(κ, p::Bool, s::system, r::Float64=s.r_max/2, algo=Vern9()) (f1, f2) = optimized_dirac_potentials(p, s) - prob_left = ODEProblem(dirac!, [0.0, 1.0], (0, r)) + prob_left = ODEProblem(dirac!, init_BC(), (0, r)) function func(E) prob_right = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r)) u_left = solve(prob_left, algo, p=(κ, E, f1, f2), saveat=[r]) From 0af47ef1613790711d989fddd9ce72552104941a Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Wed, 16 Apr 2025 16:08:32 -0400 Subject: [PATCH 6/6] Reorganize shooting method (not optional anymore) --- nucleons.jl | 54 ++++++++++++++++++++++++++--------------------------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/nucleons.jl b/nucleons.jl index 3c3f81d..393637c 100644 --- a/nucleons.jl +++ b/nucleons.jl @@ -56,37 +56,37 @@ init_BC() = [1.0, 1.0] # TODO: Why not [0.0, 1.0]? "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, - 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::Bool, E, s::system; shooting=true, normalize=true, algo=Vern9()) +function solveNucleonWf(κ, p::Bool, E, s::system; normalize=true, algo=Vern9()) (f1, f2) = optimized_dirac_potentials(p, s) - if shooting - @assert s.divs % 2 == 0 "divs must be an even number when shooting=true" - prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, s.r_max / 2)) - sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) - wf_right = reverse(hcat(sol.u...); dims=2) - next_r_max = s.r_max / 2 # for the next step - else - next_r_max = s.r_max - end - - prob = ODEProblem(dirac!, init_BC(), (0, next_r_max)) - sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=Δr(s)) - wf = hcat(sol.u...) + # partitioning + mid_idx = s.divs ÷ 2 + r_mid = rs(s)[mid_idx] + left_r = rs(s)[1:mid_idx] + right_r = rs(s)[mid_idx:end] - if shooting # join two segments - u1 = wf[:, end] - u2 = wf_right[:, 1] - if norm(u2) < 1e-10 - @warn "Right partition too small to rescale, setting to zero" - wf_right .= 0.0 - else - proj = only(u1' * u2) / norm(u2)^2 - wf_right .*= proj - end - wf = hcat(wf[:, 1:(end - 1)], wf_right) + # left partition + prob = ODEProblem(dirac!, init_BC(), (0, r_mid)) + sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=left_r) + wf_left = hcat(sol.u...) + + # right partition + prob = ODEProblem(dirac!, asymp_BC(κ, p, E, s.r_max), (s.r_max, r_mid)) + sol = solve(prob, algo, p=(κ, E, f1, f2), saveat=right_r) + wf_right = reverse(hcat(sol.u...); dims=2) + + # join two segments + u1 = wf_left[:, end] + u2 = wf_right[:, 1] + if norm(u2) < 1e-10 + @warn "Right partition too small to rescale, setting to zero" + wf_right .= 0.0 + else + proj = only(u1' * u2) / norm(u2)^2 + wf_right .*= proj end + wf = hcat(wf_left[:, 1:(end - 1)], wf_right) if normalize wf ./= norm(wf) * sqrt(Δr(s)) # integration by Reimann sum @@ -169,7 +169,7 @@ function calculateNucleonDensity(p::Bool, s::system)::Tuple{Vector{Float64}, Vec ρr2 = zeros(2, s.divs + 1) # ρ×r² values 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 ρr2 += (occ / (4 * pi)) * wf2 # 2j+1 factor is accounted in the occupancy number end