From db9dd38803ffcab4248038cf1d8cbb3afacacad2 Mon Sep 17 00:00:00 2001 From: arturgower Date: Mon, 25 Sep 2023 22:27:05 +0100 Subject: [PATCH] adjusting pair-correlation to match volume fraction --- src/ParticleCorrelations.jl | 7 +- src/inverse-problem.jl | 127 ++++++++++++++++++++ src/pair-correlation.jl | 22 +++- test/particulate-from-structure.jl | 162 ++++++++++++------------- test/particulate-hyperuniform.jl | 182 +++++++++++++++++++++++++++++ 5 files changed, 415 insertions(+), 85 deletions(-) create mode 100644 src/inverse-problem.jl create mode 100644 test/particulate-hyperuniform.jl diff --git a/src/ParticleCorrelations.jl b/src/ParticleCorrelations.jl index 35de57c..fe8ceda 100644 --- a/src/ParticleCorrelations.jl +++ b/src/ParticleCorrelations.jl @@ -1,8 +1,9 @@ module ParticleCorrelations -export Specie, Species, number_density, volume_fraction, exclusion_distance +export Specie, Species, volume_fraction, exclusion_distance export HardMedium -export periodic_particles +export number_density, translate_pair_correlation +export periodic_particles, trapezoidal_scheme, optimise_particulate ## Pair correlation export DiscretePairCorrelation, DiscreteStructureFactor @@ -32,5 +33,7 @@ include("particles.jl") include("pair-correlation.jl") include("structure-factor.jl") include("numerical.jl") +include("inverse-problem.jl") + end diff --git a/src/inverse-problem.jl b/src/inverse-problem.jl new file mode 100644 index 0000000..37fc88f --- /dev/null +++ b/src/inverse-problem.jl @@ -0,0 +1,127 @@ +function optimise_particulate(target_S::DiscreteStructureFactor{Dim}, specie::Specie{Dim}; + correlation_length::Float64 = 4.0, + numberofparticles::Int = 200, + optimoptions::Optim.Options{Float64} = Optim.Options(f_tol = 1e-5, iterations = 10) + ) where Dim + + # Define two regions R1 and R2 to place the particles + cell_number = (numberofparticles)^(1/Dim) |> round + + cell_volume = 1 / number_density(specie); + cell_length = cell_volume ^ (1/Dim) + box_length = cell_length * cell_number + + if Dim != 2 + error("the box regions only written for 2D, where as Dim = $Dim") + end + + reg1 = Box([ + [-box_length/2, -box_length/2], + [box_length/2, box_length/2] + ]); + + cell_number_boundary = ceil(correlation_length/(2 * cell_length)) + boundary_length = cell_number_boundary * cell_length + + reg2 = Box([ + [-box_length/2 - boundary_length, -box_length/2 - boundary_length], + [box_length/2 + boundary_length, box_length/2 + boundary_length] + ]); + + particles = periodic_particles(reg2,specie; random_perturbation = true) + particles1 = filter(p -> p ⊆ reg1, particles) + + S = structure_factor(particles, target_S.k; inner_box = reg1) + + points = origin.(particles) + x0 = Vector{Float64}(vcat(points...)) + + function objective(x) + points = Iterators.partition(x,Dim) |> collect + S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) + + return sum(abs2.(S.S - target_S.S)) + end + + ks = abs.(target_S.k) + + function objective_g!(G, x) + points = Iterators.partition(x,Dim) |> collect + points1 = filter(p -> p ∈ reg1, points) + J1 = length(points1) + + S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) + dS = S.S - target_S.S # add quadradture weight w here if wanted + + G[:] = vcat(map(points) do p + sum( + if p1 == p + zeros(Float64, Dim) + else + Rji = p - p1 + nRji = norm(Rji) + + Gji = (p ∈ reg1 ? 2 : 1) * sum(diffbesselj.(0, ks .* nRji) .* dS) / J1 + 2 * Gji * Rji / nRji + end + for p1 in points1) + end...) + end + + res = optimize(objective, objective_g!, x0, LBFGS(), + optimoptions + ) + + println("The result of the global step was:") + show(res) + + x0 = res.minimizer + minf = res.minimum + A = 10 * minf + 0.1 + + function penalise(x) + points = Iterators.partition(x,Dim) |> collect + α = 1 / (2*outer_radius(specie))^2 + + return A * sum( + sum( + begin + R12 = sum(abs2.(p1 - p2)) * α + if 0 < R12 < 1.0 + exp(- 4 * R12) + else 0.0 + end + end + for p2 in points) + for p1 in points) + end + + function penalise_g!(G, x) + points = Iterators.partition(x,Dim) |> collect + α = 1 / (2*outer_radius(specie))^2 + + G[:] = - (A * 16 * α) .* vcat( + map(points) do p2 + sum( + begin + R21 = sum(abs2.(p2 - p1)) .* α + if 0 < R21 < 1.0 + exp(- 4 * R21) .* (p2 - p1) + else zeros(Float64, Dim) + end + end + for p1 in points) + end...) + end + + # refine the optimisation to not allow particles to overlap + f(x) = objective(x) + penalise(x) + g(G, x) = objective_g!(G, x) + penalise_g!(G, x) + + res = optimize(f, g, x0, LBFGS(), + optimoptions + ) + + # x0 = res.minimizer + res +end \ No newline at end of file diff --git a/src/pair-correlation.jl b/src/pair-correlation.jl index d92a26d..1d73dcd 100644 --- a/src/pair-correlation.jl +++ b/src/pair-correlation.jl @@ -132,9 +132,27 @@ function number_density(discrete_pair::DiscretePairCorrelation{Dim}) where Dim sumg = sum(discrete_pair.g .* discrete_pair.r .^ (Dim - 1) .* σs) - number_density_predict = 1 / (2*(Dim - 1) * π * (discrete_pair.r[end] ^ Dim / Dim - sumg)) + term = discrete_pair.r[end] ^ Dim / Dim - sumg; -end + number_density_predict = 1 / (2*(Dim - 1) * π * term) +end + +function translate_pair_correlation(discrete_pair::DiscretePairCorrelation{Dim}, number_density::Float64) where Dim + + R = discrete_pair.r[end] + σs = trapezoidal_scheme(discrete_pair.r) + + dp = exp.(- 6rs ./ R ) + + sumg = sum(discrete_pair.g .* discrete_pair.r .^ (Dim - 1) .* σs) + sumdp = sum(dp .* discrete_pair.r .^ (Dim - 1) .* σs) + + α = (R^Dim / Dim - (1 / (2*(Dim - 1) * pi * number_density)) - sumg) / sumdp + + scaled_discrete_pair = DiscretePairCorrelation(Dim, discrete_pair.r, discrete_pair.g + α .* dp; number_density = number_density) + + return scaled_discrete_pair +end function DiscretePairCorrelation(Dim::Int,r::AbstractVector, g::AbstractVector; number_density::AbstractFloat = -1.0, diff --git a/test/particulate-from-structure.jl b/test/particulate-from-structure.jl index e32e392..24f8324 100644 --- a/test/particulate-from-structure.jl +++ b/test/particulate-from-structure.jl @@ -8,14 +8,14 @@ using Plots # let us first aim to recover a MonteCarlo pair-correlation from a crystalline arrangement # choose the spatial dimension -dim = 2 +Dim = 2 # choose the medium for the particles. Currently only one type -medium = HardMedium{dim}() +medium = HardMedium{Dim}() # choose the shapes of the particles -pairtype = MonteCarloPairCorrelation(dim; +pairtype = MonteCarloPairCorrelation(Dim; rtol = 1e-3, maxlength = 100, iterations = 10, @@ -23,14 +23,14 @@ pairtype = MonteCarloPairCorrelation(dim; ) # choose the medium for the particles. Currently only one type -medium = HardMedium{dim}() +medium = HardMedium{Dim}() # Choose the species, which represents a collection of one type of particle radius = 0.5 specie = Specie( medium, - Sphere(dim, radius), + Sphere(Dim, radius), volume_fraction = 0.15, separation_ratio = 1.0 # minimal distance from this particle = r * (separation_ratio - 1.0) ); @@ -44,87 +44,63 @@ pair = pair_correlation(specie, pairtype, rs) # If you have the Plots package plot(pair.r, pair.g) -ks = 1.0:0.5:20.0 +ks = 1.0:0.5:50.0 +ks = 1.0:0.3:40.0 +ks2 = 1.0:0.2:160.0 sfactor = structure_factor(pair, ks) +sfactor2 = structure_factor(pair, ks2) plot(sfactor.k, sfactor.S) +plot!(sfactor2.k, sfactor2.S) using Optim +# target_S = sfactor +# Dim = 2 +# numberofparticles = 200 +# correlation_length = 4 -target_S = sfactor -Dim = 2 -numberofparticles = 200 -correlation_length = 4 - -function optimise_particulate(target_S::DiscreteStructureFactor{Dim}, specie::Specie{Dim}; - correlation_length::Float64 = 4.0, - numberofparticles::Int = 200) where Dim - - # Define two regions R1 and R2 to place the particles - cell_number = (numberofparticles)^(1/Dim) |> round - - cell_volume = 1 / number_density(specie); - cell_length = cell_volume ^ (1/Dim) - box_length = cell_length * cell_number - - reg1 = Box([ - [-box_length/2, -box_length/2], - [box_length/2, box_length/2] - ]); - - cell_number_boundary = ceil(correlation_length/(2 * cell_length)) - boundary_length = cell_number_boundary * cell_length - - reg2 = Box([ - [-box_length/2 - boundary_length, -box_length/2 - boundary_length], - [box_length/2 + boundary_length, box_length/2 + boundary_length] - ]); - - particles = periodic_particles(reg2,specie; random_perturbation = true) - particles1 = filter(p -> p ⊆ reg1, particles) - - S = structure_factor(particles,ks; inner_box = reg1) - - points = origin.(particles) - x0 = vcat(points...) - - function objective(x) - points = Iterators.partition(x,Dim) |> collect - S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) - - return sum(abs2.(S.S - target_S.S)) - end - - ks = abs.(target_S.k) - - function g!(G, x) - points = Iterators.partition(x,Dim) |> collect - points1 = filter(p -> p ∈ reg1, points) - J1 = length(points1) - - S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) - dS = S.S - target_S.S # add quadradture weight w here if wanted - - G[:] = vcat(map(points) do p - sum( - if p1 == p - zeros(Float64, Dim) - else - Rji = p - p1 - nRji = norm(Rji) - - Gji = (p ∈ reg1 ? 2 : 1) * sum(diffbesselj.(0, ks .* nRji) .* dS) / J1 - 2 * Gji * Rji / nRji - end - for p1 in points1) - end...) - end - - optimize(objective, g!, x0, LBFGS()) -end +res = optimise_particulate(sfactor, specie; + numberofparticles = 300, + optimoptions = Optim.Options(g_tol = 0.0, f_tol = 1e-5, iterations = 20) +) + +show(res) + +x = res.minimizer +points = Iterators.partition(x,Dim) |> collect +x1 = [p[1] for p in points] +x2 = [p[2] for p in points] + +scatter(x1,x2) +# scatter(x1,x2, xlims =(-30.0,30.0), ylims = (-30.0,30.0)) +outer_box = Box(points); +inner_box = Box(outer_box.origin, outer_box.dimensions .- 8.0); +result_S = DiscreteStructureFactor(points, ks; + inner_box = inner_box +) +norm(result_S.S - sfactor.S) + +plot(result_S.k, result_S.S) +plot!(sfactor.k, sfactor.S, linestyle = :dash) + +result_S = DiscreteStructureFactor(points, ks2; + inner_box = inner_box +) + +plot(result_S.k, result_S.S) +plot!(sfactor2.k, sfactor2.S, linestyle = :dash) + + +result_pair = DiscretePairCorrelation(points, pair.r) +plot(pair.r, pair.g) +plot!(result_pair.r, result_pair.g) + + +points = origin.(particles) +x = vcat(points...) points = Iterators.partition(x,Dim) |> collect points1 = filter(p -> p ∈ reg1, points) @@ -133,6 +109,18 @@ J1 = length(points1) S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) dS = S.S - target_S.S # add quadradture weight w here if wanted + +function objective(x) + points = Iterators.partition(x,Dim) |> collect + S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) + + return sum(abs2.(S.S - target_S.S)) +end + +ks = abs.(target_S.k) + + + G = vcat(map(points) do p sum( if p1 == p @@ -153,22 +141,34 @@ DfDrs = Iterators.partition(G,Dim) |> collect plot(particles) plot!(reg1) -x = [p[1] for p in points]; -y = [p[2] for p in points]; +x1 = [p[1] for p in points]; +x2 = [p[2] for p in points]; u = [p[1] for p in DfDrs]; v = [p[2] for p in DfDrs]; scale = 6 -quiver!(x,y,quiver=(scale .* u,scale .* v)) +quiver!(x1,x2,quiver=(scale .* u,scale .* v)) plot!() +i_max = 40 +xs = map(1:i_max) do i + x - G .* 2i/i_max +end + +fs = objective.(xs) + +plot(fs) + +## An arteficial pair-correlation like hyper-uniform disorder. +a = 1.0 +rs = 0.2:0.4:8.0 + g = 1.0 .+ exp.(-(rs .- 4a).^2) dpair = DiscretePairCorrelation(2, rs, g) number_density(dpair) - plot(rs,g) diff --git a/test/particulate-hyperuniform.jl b/test/particulate-hyperuniform.jl new file mode 100644 index 0000000..1f3b0a4 --- /dev/null +++ b/test/particulate-hyperuniform.jl @@ -0,0 +1,182 @@ +# Here we calculate a particle configuration from a structure factor + +using ParticleCorrelations +using Test, Statistics, LinearAlgebra + +using Plots + +# let us first aim to recover a MonteCarlo pair-correlation from a crystalline arrangement + +# choose the spatial dimension +Dim = 2 + +# choose the medium for the particles. Currently only one type +medium = HardMedium{Dim}() + + +# choose the medium for the particles. Currently only one type +medium = HardMedium{Dim}() + +# Choose the species, which represents a collection of one type of particle +radius = 0.5 + +specie = Specie( + medium, + Sphere(Dim, radius), + volume_fraction = 0.15, + separation_ratio = 1.0 # minimal distance from this particle = r * (separation_ratio - 1.0) +); + +specie.particle.medium + +rs = 0.2:0.4:8.0 + +## An artificial pair-correlation like hyper-uniform disorder. +a = 1.0 +R = 8.0; +rs = (2a):0.01:R + +g = 1.0 .+ sin.(-2 .* (rs .- 3.6a) .^ 2) ./ (-2 .* (rs .- 3.6a).^ 2) +# g = 1.0 .+ cos.(-2 .* (rs .- 2a)) .* exp.(- rs) +# g = 1.0 .+ exp.(- 6rs ./ R ) + +dpair = DiscretePairCorrelation(2, rs, g) + +plot(dpair.r, dpair.g) + +number_density(dpair) + + +volfrac = 0.15; +numdensity = volfrac / volume(specie) + +dpair = translate_pair_correlation(dpair, numdensity) +plot!(dpair.r, dpair.g) +plot!([dpair.r[1],dpair.r[end]], [1.0,1.0]) + +number_density(dpair) +# pair = pair_correlation(Dim, r, g) + + +# If you have the Plots package +plot(pair.r, pair.g) + +ks = 1.0:0.5:50.0 +ks = 1.0:0.3:40.0 +ks2 = 1.0:0.2:160.0 +sfactor = structure_factor(pair, ks) +sfactor2 = structure_factor(pair, ks2) + +plot(sfactor.k, sfactor.S) +plot!(sfactor2.k, sfactor2.S) + +using Optim + +# target_S = sfactor +# Dim = 2 +# numberofparticles = 200 +# correlation_length = 4 + +res = optimise_particulate(sfactor, specie; + numberofparticles = 300, + optimoptions = Optim.Options(g_tol = 0.0, f_tol = 1e-5, iterations = 20) +) + +show(res) + +x = res.minimizer +points = Iterators.partition(x,Dim) |> collect +x1 = [p[1] for p in points] +x2 = [p[2] for p in points] + +scatter(x1,x2) +# scatter(x1,x2, xlims =(-30.0,30.0), ylims = (-30.0,30.0)) + +outer_box = Box(points); +inner_box = Box(outer_box.origin, outer_box.dimensions .- 8.0); + +result_S = DiscreteStructureFactor(points, ks; + inner_box = inner_box +) +norm(result_S.S - sfactor.S) + +plot(result_S.k, result_S.S) +plot!(sfactor.k, sfactor.S, linestyle = :dash) + +result_S = DiscreteStructureFactor(points, ks2; + inner_box = inner_box +) + +plot(result_S.k, result_S.S) +plot!(sfactor2.k, sfactor2.S, linestyle = :dash) + + +result_pair = DiscretePairCorrelation(points, pair.r) +plot(pair.r, pair.g) +plot!(result_pair.r, result_pair.g) + + +points = origin.(particles) +x = vcat(points...) + +points = Iterators.partition(x,Dim) |> collect +points1 = filter(p -> p ∈ reg1, points) +J1 = length(points1) + +S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) +dS = S.S - target_S.S # add quadradture weight w here if wanted + + +function objective(x) + points = Iterators.partition(x,Dim) |> collect + S = DiscreteStructureFactor(points, target_S.k; inner_box = reg1) + + return sum(abs2.(S.S - target_S.S)) +end + +ks = abs.(target_S.k) + + + +G = vcat(map(points) do p + sum( + if p1 == p + zeros(Float64, Dim) + else + Rji = p - p1 + nRji = norm(Rji) + + Gji = (p ∈ reg1 ? 2 : 1) * sum(diffbesselj.(0, ks .* nRji) .* dS) / J1 + 2 * Gji * Rji / nRji + end + for p1 in points1) +end...) + + +DfDrs = Iterators.partition(G,Dim) |> collect + +plot(particles) +plot!(reg1) + +x1 = [p[1] for p in points]; +x2 = [p[2] for p in points]; + +u = [p[1] for p in DfDrs]; +v = [p[2] for p in DfDrs]; + +scale = 6 +quiver!(x1,x2,quiver=(scale .* u,scale .* v)) +plot!() + +i_max = 40 +xs = map(1:i_max) do i + x - G .* 2i/i_max +end + +fs = objective.(xs) + +plot(fs) + + +plot(rs,g) +