fix instability in high d, bump minor version

Project.toml

@@ -1,24 +1,26 @@
 name = "OrnsteinZernike"
 uuid = "51399321-c693-48c0-994d-05e3df0250f2"
 authors = ["Ilian Pihlajamaa <[email protected]>"]
-version = "0.1.3"
+version = "0.1.4"
 
 [deps]
+Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+FunctionZeros = "b21f74c0-b399-568f-9643-d20f4fa2c814"
 Hankel = "74863788-d124-456e-a676-9b76578dd39e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+Dierckx = "0.5"
 FFTW = "1"
+ForwardDiff = "0.10"
 Hankel = "0.5"
+Roots = "2.0"
 SpecialFunctions = "2"
 StaticArrays = "1"
 julia = "1.7"
-Roots = "2.0"
-ForwardDiff = "0.10"
-Dierckx = "0.5"

src/Closures.jl

@@ -30,9 +30,13 @@ closure = PercusYevick()
 struct PercusYevick <: Closure end
 
 function closure_cmulr_from_gammamulr(::PercusYevick, r::Number, mayer_f::T, Γmulr::T, βuLR) where T
-    return  @. r*(mayer_f+1)*exp(βuLR)*(1+Γmulr/r-βuLR) - Γmulr - r#mayer_f*(r + Γmulr)
+    return  @. (mayer_f+1)*exp(βuLR)*(r+Γmulr-βuLR*r) - Γmulr - r#mayer_f*(r + Γmulr)
 end
 
+# function closure_c_from_gamma(::PercusYevick, r::Number, mayer_f::T, Γmulr::T, βuLR) where T
+#     return  @. mayer_f*(r + Γmulr)
+# end
+
 function bridge_function(::PercusYevick, _, _, γ)
     B  = @. log1p(γ) - γ
     return B

src/FourierTransforms.jl

@@ -15,6 +15,16 @@ struct My1DPlan{T2,T3,T4, T5, T6}
     M::T4
 end
 
+struct MyNDPlan{dims, T, T2}
+    r::Vector{T}
+    k::Vector{T}
+    M::Int
+    plan_forward::Matrix{T2}
+    plan_backward::Matrix{T2}
+end
+
+
+
 function get_fourier_plan(::SimpleLiquid{1, species, T1, T2, P}, method, F) where {species, T1, T2, P}
     M = length(F)
     dr = method.dr
@@ -65,6 +75,27 @@ function find_fourier_plan_nd(::SimpleLiquid{ndims, species, T1, T2, P}, F::Vect
     M = length(F)
     plan = QDHT(ndims/2-1, 1, M*dr, M)
     return plan
+    # M = length(F)
+    # plan_forward = zeros(eltype(eltype(F)), M, M)
+    # plan_backward = zeros(eltype(eltype(F)), M, M)
+    # println("hi")
+    # p = ndims/2 - 1
+    # λ_i = besselj_zero.(p, 1:(M+1))
+    # Rmax = dr * M
+    # Kmax = λ_i[end]/Rmax
+    # k = λ_i[1:end-1]/Rmax
+    # r = λ_i[1:end-1]/Kmax
+
+    # for i = 1:M
+    #     for j = 1:M
+    #         plan_forward[j, i] = 2 * (2π)^(p+1) / Kmax^2 *  besselj(p, k[j]*r[i])/besselj(p+1, Kmax*r[i])^2
+    #         plan_backward[j, i] = 2 / ( Rmax^2 * (2π)^(p+1) ) * besselj(p, k[i]*r[j])/besselj(p+1, Rmax*k[i])^2
+    #     end
+    # end
+    # Tr = eltype(r)
+    # Tq = eltype(plan_backward)
+    # return MyNDPlan{ndims, Tr, Tq}(r, k, M, plan_forward, plan_backward)
+
 end
 
 function inverse_radial_fourier_transform_3d(F̂, r, k)
@@ -181,6 +212,17 @@ function fourier!(F̂::AbstractVector{T}, F::AbstractVector{T}, Q::Hankel.QDHT{p
     @. F = F / r^(halfdims-2)
 
     @. F̂ = F̂ * (2π)^halfdims / k^(halfdims-2)
+
+    # p = dims/2 - 1
+    # k = Q.k
+    # r = Q.r
+
+    # @. F = F * r^(p-1)
+    # mul!(F̂, Q.plan_forward, F)
+    # @. F = F / r^(p-1)
+
+    # @. F̂ = F̂ / k^(p-1)
+
 end
 
 """
@@ -199,6 +241,15 @@ function inverse_fourier!(F::AbstractVector{T}, F̂::AbstractVector{T}, Q::Hanke
     @. F̂ = F̂ / k^(halfdims-2)
     @. F ./= Q.scaleRK^2
     @. F = F * (2π)^(-halfdims) / r^(halfdims-2)
+
+    # p = dims/2 - 1
+    # k = Q.k
+    # r = Q.r
+
+    # @. F̂ = F̂ * k^(p-1)
+    # mul!(F, Q.plan_backward, F̂)
+    # @. F̂ = F̂ / k^(p-1)
+    # @. F = F / r^(p-1)
 end
 
 

src/OrnsteinZernike.jl

@@ -5,6 +5,8 @@ A generic solver package for Ornstein-Zernike equations from liquid state theory
 """
 module OrnsteinZernike
 using FFTW, StaticArrays, LinearAlgebra, Hankel, SpecialFunctions, Dierckx 
+using Bessels: besselj
+using FunctionZeros
 using Roots: find_zero
 
 export solve

src/Solvers/FourierIteration.jl

@@ -57,7 +57,7 @@ function solve(system::SimpleLiquid{dims, species, T1, T2, P}, closure::Closure,
     end
     if method.verbose 
         print("Converged after $iteration iterations, ")
-        println("the error is $(round(err, digits=ceil(Int, 1-log10(tolerance)))).")
+        println("the error is $(round(err, digits=ceil(Int, 1-log10(err)))).")
     end
     c = C ./ r
     g = find_g_from_c_and_Γ(c, Γ_new, r)

src/Solvers/NgIteration.jl

@@ -81,7 +81,7 @@ function solve(system::SimpleLiquid{dims, 1, T1, T2, P}, closure::Closure, metho
 
     if method.verbose 
         print("Converged after $iteration iterations, ")
-        println("the error is $(round(err, digits=ceil(Int, 1-log10(tolerance)))).")
+        println("the error is $(round(err, digits=ceil(Int, 1-log10(err)))).")
     end
     c = C ./ r
     g = find_g_from_c_and_Γ(c, Γ_new, r)
@@ -189,7 +189,7 @@ function solve(system::SimpleLiquid{dims, species, T1, T2, P}, closure::Closure,
 
     if method.verbose 
         print("Converged after $iteration iterations, ")
-        println("the error is $(round(err, digits=ceil(Int, 1-log10(tolerance)))).")
+        println("the error is $(round(err, digits=ceil(Int, 1-log10(err)))).")
     end
     c = C ./ r
     g = find_g_from_c_and_Γ(c, Γ_new, r)