diff options
Diffstat (limited to 'src/anyeq.jl')
| -rw-r--r-- | src/anyeq.jl | 1485 |
1 files changed, 1196 insertions, 289 deletions
diff --git a/src/anyeq.jl b/src/anyeq.jl index 8459cb0..fee77d8 100644 --- a/src/anyeq.jl +++ b/src/anyeq.jl @@ -1,4 +1,4 @@ -## Copyright 2021 Ian Jauslin +## Copyright 2021-2023 Ian Jauslin ## ## Licensed under the Apache License, Version 2.0 (the "License"); ## you may not use this file except in compliance with the License. @@ -28,24 +28,26 @@ end # compute energy for a given rho -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) @@ -53,72 +55,60 @@ function anyeq_energy(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx end # compute energy as a function of rho -function anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - - # save result from each task - es=Array{Float64,1}(undef,length(rhos)) - err=Array{Float64,1}(undef,length(rhos)) + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) - ## spawn workers - # number of workers - nw=nworkers() - # split jobs among workers - work=Array{Array{Int64,1},1}(undef,nw) - # init empty arrays - for p in 1:nw - work[p]=zeros(0) - end - for j in 1:length(rhos) - append!(work[(j-1)%nw+1],j) - end + # spawn workers + work=spawn_workers(length(rhos)) - count=0 - # for each worker - @sync for p in 1:nw - # for each task - @async for j in work[p] - count=count+1 - if count>=nw - progress(count,length(rhos),10000) - end - # run the task - (u,es[j],err[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) - end - end + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) for j in 1:length(rhos) - @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j]) + @printf("% .15e % .15e % .15e\n",rhos[j],es[j],errs[j]) end end # compute energy as a function of rho # initialize with previous rho -function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy_rho_init_prevrho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - u0s=anyeq_init_medeq([rhos[1]],N,J,k,a0,v,maxiter,tolerance) - u=u0s[1] + u=anyeq_init_medeq([rhos[1]],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) for j in 1:length(rhos) progress(j,length(rhos),10000) @@ -134,20 +124,26 @@ function anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,ap end # compute energy as a function of rho # initialize with next rho -function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_energy_rho_init_nextrho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - u0s=anyeq_init_medeq([rhos[length(rhos)]],N,J,k,a0,v,maxiter,tolerance) - u=u0s[1] + u=anyeq_init_medeq([rhos[length(rhos)]],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) for j in 1:length(rhos) progress(j,length(rhos),10000) @@ -163,23 +159,26 @@ function anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,ap end # compute u(k) -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_uk( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # init vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) @@ -192,58 +191,60 @@ function anyeq_uk(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,approx,sav end # compute u(x) -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_ux(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile) +function anyeq_ux( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Anyeq_approx, + savefile::String +) # init vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) for i in 1:nx x=xmin+(xmax-xmin)*i/nx - ux=0. - for zeta in 0:J-1 - for j in 1:N - ux+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*u[zeta*N+j]*sin(k[zeta*N+j]*x) - end - end - @printf("% .15e % .15e % .15e\n",x,real(ux),imag(ux)) + ux=anyeq_u_x(x,u,k,N,J,taus,weights) + @printf("% .15e % .15e % .15e\n",x,ux,) end end # compute condensate fraction for a given rho -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) +function anyeq_condensate_fraction( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) - end + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) @@ -254,161 +255,571 @@ function anyeq_condensate_fraction(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tol end # condensate fraction as a function of rho -function anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) - ## spawn workers - # number of workers - nw=nworkers() - # split jobs among workers - work=Array{Array{Int64,1},1}(undef,nw) - # init empty arrays - for p in 1:nw - work[p]=zeros(0) +function anyeq_condensate_fraction_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) + # spawn workers + work=spawn_workers(length(rhos)) + + # initialize vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) + + # compute eta + etas=Array{Float64,1}(undef,length(rhos)) + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + etas[j]=remotecall_fetch(anyeq_eta,workers()[p],us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + end end + for j in 1:length(rhos) - append!(work[(j-1)%nw+1],j) + @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j]) end +end +# compute the momentum distribution for a given rho +function anyeq_momentum_distribution( + kmin::Float64, + kmax::Float64, + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + # compute M + if windowL==Inf + # use discrete approximation of delta function + M=anyeq_momentum_discrete_delta(kmin,kmax,u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + M=anyeq_momentum_window(kmin,kmax,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) end + # order k's in increasing order + for zeta in 0:J-1 + for j in 1:N + q=k[(J-1-zeta)*N+j] + # drop if not in k interval + if q<kmin || q>kmax + continue + end + @printf("% .15e % .15e\n",q,M[(J-1-zeta)*N+j]) + end + end +end + +# 2 point correlation function +function anyeq_2pt_correlation( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + # compute initial guess from medeq - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) + @printf("% .15e % .15e\n",x,C2) + end +end + +# maximum of 2 point correlation function +function anyeq_2pt_correlation_max( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + dx::Float64, + x0::Float64, # initial guess is x0/rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # initial guess + x0=1/rho^(1/3) + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (x,f)=anyeq_2pt_max(u,P,N,J,x0/rho^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rho,weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + if(x==Inf) + @printf(stderr,"max search failed for rho=%e\n",rho) + else + @printf("% .15e % .15e\n",x,f) + end +end + +# maximum of 2 point correlation function as a function of rho +function anyeq_2pt_correlation_max_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + dx::Float64, + x0::Float64, # initial guess is x0/rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # save result from each task + xs=Array{Float64,1}(undef,length(rhos)) + fs=Array{Float64,1}(undef,length(rhos)) + + # spawn workers + work=spawn_workers(length(rhos)) # compute u - us=Array{Array{Float64,1}}(undef,length(rhos)) - errs=Array{Float64,1}(undef,length(rhos)) + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) + count=0 # for each worker - @sync for p in 1:nw + @sync for p in 1:length(work) # for each task @async for j in work[p] count=count+1 - if count>=nw - progress(count,length(rhos),10000) + if count>=length(work) + progress(count,length(rhos),10000) end # run the task - (us[j],E,errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (xs[j],fs[j])=remotecall_fetch(anyeq_2pt_max,workers()[p],us[j],P,N,J,x0/rhos[j]^(1/3),dx,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) end end - # compute eta - etas=Array{Float64}(undef,length(rhos)) + for j in 1:length(rhos) + if(xs[j]==Inf) + @printf(stderr,"max search failed for rho=%e\n",rhos[j]) + else + @printf("% .15e % .15e % .15e\n",rhos[j],xs[j],fs[j]) + end + end +end + +# Correlation function in Fourier space +function anyeq_2pt_correlation_fourier( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + kmin::Float64, + kmax::Float64, + nk::Int64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + # save result from each task + C2s=Array{Float64,1}(undef,nk) + + # spawn workers + work=spawn_workers(nk) + count=0 # for each worker - @sync for p in 1:nw + @sync for p in 1:length(work) # for each task @async for j in work[p] count=count+1 - if count>=nw - progress(count,length(rhos),10000) + if count>=length(work) + progress(count,nk,10000) end # run the task - etas[j]=anyeq_eta(us[j],P,N,J,rhos[j],weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + C2s[j]=remotecall_fetch(anyeq_2pt_fourier,workers()[p],kmin+(kmax-kmin)*j/nk,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) end end - for j in 1:length(rhos) - @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[j]) + for j in 1:nk + @printf("% .15e % .15e\n",kmin+(kmax-kmin)*j/nk,C2s[j]) end end -# compute the momentum distribution for a given rho -# use minlrho, nlrho to incrementally compute the solution to medeq (to initialize the Newton algorithm) -function anyeq_momentum_distribution(rho,minlrho,nlrho,taus,P,N,J,a0,v,maxiter,tolerance,approx,savefile) - # initialize vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) +# maximum of Fourier transform of 2 point correlation function +function anyeq_2pt_correlation_fourier_max( + rho::Float64, + minlrho_init::Float64, + nlrho_init::Int64, + dk::Float64, + k0::Float64, # initial guess is k0*rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (ko,f)=anyeq_2pt_fourier_max(u,P,N,J,k0*rho^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rho,weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + if(ko==Inf) + @printf(stderr,"max search failed for rho=%e\n",rho) else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + @printf("% .15e % .15e\n",ko,f) end +end + +# maximum of fourier transform of 2 point correlation function as a function of rho +function anyeq_2pt_correlation_fourier_max_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + dk::Float64, + k0::Float64, # initial guess is k0*rho^(1/3) + maxstep::Float64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + tolerance_max::Float64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) - end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) - (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + # spawn workers + work=spawn_workers(length(rhos)) - # compute M - M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) - for zeta in 0:J-1 - for j in 1:N - # order k's in increasing order - @printf("% .15e % .15e\n",k[(J-1-zeta)*N+j],M[(J-1-zeta)*N+j]) + # save result from each task + ks=Array{Float64,1}(undef,length(rhos)) + fs=Array{Float64,1}(undef,length(rhos)) + + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + (ks[j],fs[j])=remotecall_fetch(anyeq_2pt_fourier_max,workers()[p],us[j],P,N,J,k0*rhos[j]^(1/3),dk,maxstep,maxiter,tolerance_max,windowL,rhos[j],weights,T,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + end + end + + for j in 1:length(rhos) + if(ks[j]==Inf) + @printf(stderr,"max search failed for rho=%e\n",rhos[j]) + else + @printf("% .15e % .15e % .15e\n",rhos[j],ks[j],fs[j]) end end end -# 2 point correlation function -function anyeq_2pt_correlation(minlrho,nlrho,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx,savefile) +# uncondensed 2 point correlation function +function anyeq_uncondensed_2pt_correlation( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmin::Float64, + xmax::Float64, + nx::Int64, + approx::Anyeq_approx, + savefile::String +) # init vectors - (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) - # init Abar - if savefile!="" - Abar=anyeq_read_Abar(savefile,P,N,J) - else - Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # compute u and some useful integrals + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + # compute u(xi) + uxis=Array{Float64,1}(undef,nx) + for i in 1:nx + uxis[i]=anyeq_u_x(xmin+(xmax-xmin)*i/nx,u,k,N,J,taus,weights) + end + + + # spawn workers + work=spawn_workers(nx) + gamma2s=Array{Float64,1}(undef,nx) + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,nx,10000) + end + # run the task + gamma2s[j]=remotecall_fetch(anyeq_uncondensed_2pt,workers()[p],xmin+(xmax-xmin)*j/nx,uxis[j],u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) + end end + for i in 1:nx + @printf("% .15e % .15e\n",xmin+(xmax-xmin)*i/nx,gamma2s[i]) + end +end + +# compute the compressibility +function anyeq_compressibility_rho( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + savefile::String +) + # initialize vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + # compute initial guess from medeq - rhos=Array{Float64}(undef,nlrho) - for j in 0:nlrho-1 - rhos[j+1]=(nlrho==1 ? rho : 10^(minlrho+(log10(rho)-minlrho)/(nlrho-1)*j)) + u0s=anyeq_init_medeq(rhos,minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) + + # save result from each task + es=Array{Float64,1}(undef,length(rhos)) + err=Array{Float64,1}(undef,length(rhos)) + + # spawn workers + work=spawn_workers(length(rhos)) + + # compute u + (us,es,errs)=anyeq_hatu_rho_multithread(u0s,P,N,J,rhos,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx,work) + + for j in 1:length(rhos)-2 + # compute \rho^2\partial^2(\rho e)=\partial^2_{\log\rho}(\rho e)-\partial_{\log\rho}(\rho e) as a function of rho + # \partial^2_{\log\rho}(\rho e) + p2=(rhos[j+2]*es[j+2]-2*rhos[j+1]*es[j+1]+rhos[j]*es[j])/(log(rhos[j+2])-log(rhos[j+1]))/(log(rhos[j+1])-log(rhos[j])) + # \partial_{\log\rho}(\rho e) (take average of front and back) + p11=(rhos[j+2]*es[j+2]-rhos[j+1]*es[j+1])/(log(rhos[j+2])-log(rhos[j+1])) + p12=(rhos[j+1]*es[j+1]-rhos[j]*es[j])/(log(rhos[j+1])-log(rhos[j])) + # compressibility is 1/this + @printf("% .15e % .15e\n",rhos[j+1],1/(p2-(p11+p12)/2)) end - u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - u0=u0s[nlrho] +end + +# Fourier transform 2 point correlation function test: compute by transforming anyeq_2pt +function anyeq_2pt_correlation_fourier_test( + minlrho_init::Float64, + nlrho_init::Int64, + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64, + xmax::Float64, + kmin::Float64, + kmax::Float64, + nk::Int64, + approx::Anyeq_approx, + savefile::String +) + # init vectors + (weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v,approx,savefile) + + # compute initial guess from medeq + u0=anyeq_init_medeq([rho],minlrho_init,nlrho_init,N,J,k,a0,v,maxiter,tolerance) # compute u and some useful integrals (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) - for i in 1:nx - x=xmin+(xmax-xmin)*i/nx - C2=anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G) - @printf("% .15e % .15e\n",x,C2) + # compute C2 in real space + C2=Array{Float64,1}(undef,N*J) + for i in 1:N*J + if k[i]<xmax + C2[i]=anyeq_2pt(k[i],u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK)-rho^2 + else + C2[i]=0. + end + end + # invert Fourier transform + for i in 1:nk + k0=kmin+(kmax-kmin)*i/nk + hatC2=inverse_fourier_chebyshev(C2,k0,k,taus,weights,N,J) + @printf("% .15e % .15e\n",k0,hatC2) end end # compute Abar -function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) +function anyeq_Abar_multithread( + k::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + taus::Array{Float64,1}, + T::Array{Polynomial,1}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx +) if approx.bL3==0. - return [] + # empty array + return Array{Array{Float64,5}}(undef,0) end out=Array{Array{Float64,5}}(undef,J*N) - ## spawn workers - # number of workers - nw=nworkers() - # split jobs among workers - work=Array{Array{Int64,1},1}(undef,nw) - # init empty arrays - for p in 1:nw - work[p]=zeros(0) - end - for j in 1:N*J - append!(work[(j-1)%nw+1],j) - end + # spawn workers + work=spawn_workers(N*J) count=0 # for each worker - @sync for p in 1:nw + @sync for p in 1:length(work) # for each task @async for j in work[p] count=count+1 - if count>=nw + if count>=length(work) progress(count,N*J,10000) end # run the task @@ -419,14 +830,23 @@ function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) return out end + # initialize computation -@everywhere function anyeq_init(taus,P,N,J,v) +@everywhere function anyeq_init( + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + v::Function, + approx::Anyeq_approx, + savefile::String +) # Gauss-Legendre weights weights=gausslegendre(N) # initialize vectors V,k - V=Array{Float64}(undef,J*N) - k=Array{Float64}(undef,J*N) + V=Array{Float64,1}(undef,J*N) + k=Array{Float64,1}(undef,J*N) for zeta in 0:J-1 for j in 1:N xj=weights[1][j] @@ -437,7 +857,7 @@ end end end # potential at 0 - V0=v(0) + V0=v(0.) # initialize matrix A T=chebyshev_polynomials(P) @@ -449,39 +869,62 @@ end Upsilon=Upsilonmat(k,v,weights_plus) Upsilon0=Upsilon0mat(k,v,weights_plus) - return(weights,T,k,V,V0,A,Upsilon,Upsilon0) + # init Abar + if savefile!="" + Abar=anyeq_read_Abar(savefile,P,N,J) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + return(weights,T,k,V,V0,A,Abar,Upsilon,Upsilon0) end # compute initial guess from medeq -@everywhere function anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) - us_medeq=Array{Array{Float64,1}}(undef,length(rhos)) - u0s=Array{Array{Float64,1}}(undef,length(rhos)) - +@everywhere function anyeq_init_medeq( + rhos::Array{Float64,1}, + minlrho_init::Float64, + nlrho_init::Int64, + N::Int64, + J::Int64, + k::Array{Float64,1}, + a0::Float64, + v::Function, + maxiter::Int64, + tolerance::Float64 +) weights_medeq=gausslegendre(N*J) + (u0s_medeq,es,errs)=easyeq_compute_u_prevrho(rhos,minlrho_init,nlrho_init,a0,N*J,v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) - (us_medeq[1],E,err)=easyeq_hatu(easyeq_init_u(a0,J*N,weights_medeq),J*N,rhos[1],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) - u0s[1]=easyeq_to_anyeq(us_medeq[1],weights_medeq,k,N,J) - if err>tolerance - print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[1],"\n") - end - - for j in 2:length(rhos) - (us_medeq[j],E,err)=easyeq_hatu(us_medeq[j-1],J*N,rhos[j],v,maxiter,tolerance,weights_medeq,Easyeq_approx(1.,1.)) - u0s[j]=easyeq_to_anyeq(us_medeq[j],weights_medeq,k,N,J) - - if err>tolerance + # check errs + for j in 1:length(errs) + if errs[j]>tolerance print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n") end end - return u0s + # return a single vector if there is a single rho + if length(rhos)>1 + u0s=Array{Array{Float64,1}}(undef,length(rhos)) + for j in 1:length(u0s_medeq) + u0s[j]=easyeq_to_anyeq(u0s_medeq[j],weights_medeq,k,N,J) + end + return u0s + else + return easyeq_to_anyeq(u0s_medeq,weights_medeq,k,N,J) + end end # interpolate the solution of medeq to an input for anyeq -@everywhere function easyeq_to_anyeq(u_simple,weights,k,N,J) +@everywhere function easyeq_to_anyeq( + u_simple::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + N::Int64, + J::Int64 +) # reorder u_simple, which is evaluated at (1-x_j)/(1+x_j) with x_j\in[-1,1] u_s=zeros(Float64,length(u_simple)) - k_s=Array{Float64}(undef,length(u_simple)) + k_s=Array{Float64,1}(undef,length(u_simple)) for j in 1:length(u_simple) xj=weights[1][j] k_s[length(u_simple)-j+1]=(1-xj)/(1+xj) @@ -501,7 +944,27 @@ end # compute u using chebyshev expansions -@everywhere function anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) +@everywhere function anyeq_hatu( + u0::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + a0::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx +) # init # rescale by rho (that's how u is defined) u=rho*u0 @@ -512,7 +975,7 @@ end # run Newton algorithm for i in 1:maxiter-1 (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) - new=u-inv(anyeq_DXi(u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y) + new=u-inv(anyeq_DXi(u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_Xi(u,X,Y) error=norm(new-u)/norm(u) if(error<tolerance) @@ -528,8 +991,60 @@ end end +# compute u for various rho +function anyeq_hatu_rho_multithread( + u0s::Array{Array{Float64,1},1}, + P::Int64, + N::Int64, + J::Int64, + rhos::Array{Float64,1}, + a0::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + v::Function, + maxiter::Int64, + tolerance::Float64, + approx::Anyeq_approx, + work::Array{Array{Int64,1},1} +) + # compute u + us=Array{Array{Float64,1}}(undef,length(rhos)) + es=Array{Float64,1}(undef,length(rhos)) + errs=Array{Float64,1}(undef,length(rhos)) + count=0 + # for each worker + @sync for p in 1:length(work) + # for each task + @async for j in work[p] + count=count+1 + if count>=length(work) + progress(count,length(rhos),10000) + end + # run the task + (us[j],es[j],errs[j])=remotecall_fetch(anyeq_hatu,workers()[p],u0s[j],P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + end + end + + return (us,es,errs) +end + + # save Abar -function anyeq_save_Abar(taus,P,N,J,v,approx) +function anyeq_save_Abar( + taus::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + v::Function, + approx::Anyeq_approx +) # initialize vectors (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) @@ -555,7 +1070,12 @@ function anyeq_save_Abar(taus,P,N,J,v,approx) end # read Abar -function anyeq_read_Abar(savefile,P,N,J) +function anyeq_read_Abar( + savefile::String, + P::Int64, + N::Int64, + J::Int64 +) # open file ff=open(savefile,"r") # read all lines @@ -623,13 +1143,39 @@ end # Xi # takes the vector of kj's and xn's as input -@everywhere function anyeq_Xi(U,X,Y) +@everywhere function anyeq_Xi( + U::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1} +) return U-(Y.+1)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2)) end # DXi # takes the vector of kj's as input -@everywhere function anyeq_DXi(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) +@everywhere function anyeq_DXi( + U::Array{Float64,1}, + rho::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) out=Array{Float64,2}(undef,N*J,N*J) for zetapp in 0:J-1 for n in 1:N @@ -720,7 +1266,23 @@ end end # compute S,E,I,J,X and Y -@everywhere function anyeq_SEIJGXY(U,rho,k,taus,v,v0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) +@everywhere function anyeq_SEIJGXY( + U::Array{Float64,1}, + rho::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + v::Array{Float64,1}, + v0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx +) # Chebyshev expansion of U FU=chebyshev(U,taus,weights,P,N,J,2) @@ -798,79 +1360,207 @@ end return(S,E,II,JJ,X,Y,sL1,sK,G) end + +# u(x) +@everywhere function anyeq_u_x( + x::Float64, + u::Array{Float64,1}, + k::Array{Float64,1}, + N::Int64, + J::Int64, + taus::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}} +) + ux=0. + for zeta in 0:J-1 + for j in 1:N + ux+=(taus[zeta+2]-taus[zeta+1])/(16*pi*x)*weights[2][j]*cos(pi*weights[1][j]/2)*(1+k[zeta*N+j])^2*k[zeta*N+j]*u[zeta*N+j]*sin(k[zeta*N+j]*x) + end + end + return real(ux) +end + # condensate fraction -@everywhere function anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) +@everywhere function anyeq_eta( + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) # compute dXi/dmu (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) dXidmu=(Y.+1)./(rho*sL1)./(2*(X.+1).^2).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1).^2 ./((X.+1).^4)./(rho*sL1).*dotdPhi((Y.+1)./(X.+1).^2) # compute eta - du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu + du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu eta=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 return eta end -# momentum distribution -@everywhere function anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) - # compute dXi/dlambda (without delta functions) - (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) - dXidlambda=-(2*pi)^3*2*u./sL1.*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2)) - - # approximation for delta function (without Kronecker deltas) - delta=Array{Float64}(undef,J*N) - for zeta in 0:J-1 - for n in 1:N - delta[zeta*N+n]=2/pi^2/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2*k[zeta*N+n]^2) - end - end - - # compute dXidu - dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)) - M=Array{Float64}(undef,J*N) - for i in 1:J*N - # du/dlambda - du=dXidu[:,i]*dXidlambda[i]*delta[i] +# correlation function +@everywhere function anyeq_2pt( + x::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) + g=(r,x)->sinc(r*x)*hann(r,windowL) + du=anyeq_dudv(g, x, u, P, N, J, rho, weights, k, taus, V, V0, A, Abar, Upsilon, Upsilon0, approx, S, E, II, JJ, X, Y, sL1, sK) + + C2=rho^2*(1-integrate_f_chebyshev(s->g(s,x),u,k,taus,weights,N,J)-integrate_f_chebyshev(s->1.,V.*du,k,taus,weights,N,J)) +end - # compute M - M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 +# uncondensed correlation function +@everywhere function anyeq_uncondensed_2pt( + xi::Float64, + uxi::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) + # compute dXi/dmu + g=Array{Float64,1}(undef,length(k)) + for i in 1:length(k) + g[i]=-2*rho*uxi*sinc(k[i]*xi)*hann(k[i],windowL) end + dXidmu=-g./sL1./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))-(Y.+1).*g./sL1./(2(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2) - return M -end + # compute gamma2 + du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*dXidmu + gamma2=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 + return gamma2 +end -# correlation function -@everywhere function anyeq_2pt(x,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK,G) +# compute the directional derivative of u with respect to v in direction g +@everywhere function anyeq_dudv( + g::Function,# should be of the form g(k,x) where x is a parameter + x::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) # initialize dV - dV=Array{Float64}(undef,J*N) + dV=Array{Float64,1}(undef,J*N) for i in 1:J*N - if x>0 - dV[i]=sin(k[i]*x)/(k[i]*x)*hann(k[i],windowL) - else - dV[i]=hann(k[i],windowL) - end + dV[i]=g(k[i],x) end - dV0=1. + dV0=g(0.,x) # compute dUpsilon # Upsilonmat does not use splines, so increase precision weights_plus=gausslegendre(N*J) - dUpsilon=Upsilonmat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus) - dUpsilon0=Upsilon0mat(k,r->sin(r*x)/(r*x)*hann(r,windowL),weights_plus) + dUpsilon=Upsilonmat(k,s->g(s,x),weights_plus) + dUpsilon0=Upsilon0mat(k,s->g(s,x),weights_plus) - du=-inv(anyeq_DXi(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(x,rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) + du=-inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK))*anyeq_dXidv(rho*u,rho,k,taus,dV,dV0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) # rescale rho du=du/rho - C2=rho^2*(1-integrate_f_chebyshev(s->1.,u.*dV+V.*du,k,taus,weights,N,J)) - - return C2 + return du end -# derivative of Xi with respect to v in the direction sin(kx)/kx -@everywhere function anyeq_dXidv(x,U,rho,k,taus,dv,dv0,A,Abar,dUpsilon,dUpsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK) +# derivative of Xi with respect to v in the specified by dUpsilon and dUpsilon0 +@everywhere function anyeq_dXidv( + U::Array{Float64,1}, + rho::Float64, + k::Array{Float64,1}, + taus::Array{Float64,1}, + dv::Array{Float64,1}, + dv0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + dUpsilon::Array{Array{Float64,1},1}, + dUpsilon0::Array{Float64,1}, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + P::Int64, + N::Int64, + J::Int64, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) # Chebyshev expansion of U FU=chebyshev(U,taus,weights,P,N,J,2) @@ -896,7 +1586,7 @@ end dJJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,dFS,Abar) end if approx.bL3!=1. - dJJ=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2 + dJJ+=approx.gL3*(1-approx.bL3)*dE*(UU/rho).^2 end end @@ -947,3 +1637,220 @@ end out=((Y.+1).*dX./(X.+1)-dY)./(2*(X.+1)).*dotPhi((Y.+1)./((X.+1).^2))+(Y.+1)./(2*(X.+1).^3).*(2*(Y.+1)./(X.+1).*dX-dY).*dotdPhi((Y.+1)./(X.+1).^2) return out end + +# maximum of 2 point correlation function +@everywhere function anyeq_2pt_max( + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + x0::Float64, + dx::Float64, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + T::Array{Polynomial,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + # compute some useful integrals + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + (x,f)=newton_maximum(y->anyeq_2pt(y,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK),x0,dx,maxiter,tolerance,maxstep) + + return(x,f) +end + + +# Fourier transform of 2pt correlation function at q +@everywhere function anyeq_2pt_fourier( + q::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx, + S::Array{Float64,1}, + E::Float64, + II::Array{Float64,1}, + JJ::Array{Float64,1}, + X::Array{Float64,1}, + Y::Array{Float64,1}, + sL1::Array{Float64,1}, + sK::Array{Float64,1} +) + # direction in which to differentiate u + weights_plus=gausslegendre(N*J) + #g=(r,x)->(r>0. ? (1.)/(2*x*r)*integrate_legendre(s->s*gaussian(s,(1.)/windowL),abs(x-r),x+r,weights_plus) : gaussian(x,(1.)/windowL)) + g=(r,x)->(r>0. ? (1.)/(2*x*r*windowL)*(gaussian(x-r,(1.)/windowL)-gaussian(x+r,(1.)/windowL)) : gaussian(x,(1.)/windowL)) + + du=anyeq_dudv(g,q,u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK) + + C2=rho^2*(-integrate_f_chebyshev(s->g(s,q),u,k,taus,weights,N,J)-integrate_f_chebyshev(s->1.,V.*du,k,taus,weights,N,J)) + + return C2 +end + +# maximum of Fourier transform of 2 point correlation function +@everywhere function anyeq_2pt_fourier_max( + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + k0::Float64, + dk::Float64, + maxstep::Float64, + maxiter::Int64, + tolerance::Float64, + windowL::Float64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + T::Array{Polynomial,1}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + # compute some useful integrals + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + + (ko,f)=newton_maximum(y->anyeq_2pt_fourier(y,u,P,N,J,windowL,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx,S,E,II,JJ,X,Y,sL1,sK),k0,dk,maxiter,tolerance,maxstep) + + return(ko,f) +end + +# momentum distribution, computed using a Gaussian window +@everywhere function anyeq_momentum_window( + kmin::Float64, + kmax::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + windowL::Float64, # L is windowL/k^2 + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + + # compute dXi/dlambda without the delta function of u(q) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + dXidlambda=(16*pi^3)*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2))./sL1 + + # compute dXidu + dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)) + + M=Array{Float64,1}(undef,N*J) + for j in 1:J*N + # the Chebyshev polynomial expansion is often not good enough to compute M away from k[i] + q=k[j] + + # drop the computation if not in k interval + if q<kmin || q>kmax + continue + end + + # delta function + delta=Array{Float64,1}(undef,J*N) + L=windowL/q^2 + for i in 1:J*N + delta[i]=(1.)/(2*k[i]*q*L)*(gaussian(k[i]-q,(1.)/L)-gaussian(k[i]+q,(1.)/L)) + end + + # du/dlambda + du=-dXidu*(dXidlambda.*delta*u[j]) + # rescale u + du=du/rho + + # compute M + M[j]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)*rho/2 + end + + return M +end + +# momentum distribution, computed using a discrete approximation of the delta function +@everywhere function anyeq_momentum_discrete_delta( + kmin::Float64, + kmax::Float64, + u::Array{Float64,1}, + P::Int64, + N::Int64, + J::Int64, + rho::Float64, + weights::Tuple{Array{Float64,1},Array{Float64,1}}, + k::Array{Float64,1}, + taus::Array{Float64,1}, + V::Array{Float64,1}, + V0::Float64, + A::Array{Array{Float64,2},1}, + Abar::Array{Array{Float64,5},1}, + Upsilon::Array{Array{Float64,1},1}, + Upsilon0::Array{Float64,1}, + approx::Anyeq_approx +) + # compute dXi/dlambda (without delta functions) + (S,E,II,JJ,X,Y,sL1,sK,G)=anyeq_SEIJGXY(rho*u,rho,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx) + dXidlambda=(2*pi)^3*2*u./sL1.*(dotPhi((Y.+1)./((X.+1).^2))./(2*(X.+1))+(Y.+1)./(2*(X.+1).^3).*dotdPhi((Y.+1)./(X.+1).^2)) + + # approximation for delta function (without Kronecker deltas) + delta=Array{Float64,1}(undef,J*N) + for zeta in 0:J-1 + for n in 1:N + delta[zeta*N+n]=8/pi/((taus[zeta+2]-taus[zeta+1])*weights[2][n]*cos(pi*weights[1][n]/2)*(1+k[zeta*N+n])^2) + end + end + + # compute dXidu + dXidu=inv(anyeq_DXi(rho*u,rho,k,taus,A,Abar,Upsilon,Upsilon0,weights,P,N,J,approx,S,E,II,JJ,X,Y,sL1,sK)) + + M=Array{Float64,1}(undef,J*N) + for i in 1:J*N + # drop the computation if not in k interval + if k[i]<kmin || k[i]>kmax + continue + end + + # du/dlambda + du=-dXidu[:,i]*dXidlambda[i]*delta[i] + + # compute M + M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 + end + + return M +end |