diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/anyeq.jl | 949 | ||||
| -rw-r--r-- | src/chebyshev.jl | 279 | ||||
| -rw-r--r-- | src/easyeq.jl | 411 | ||||
| -rw-r--r-- | src/integration.jl | 58 | ||||
| -rw-r--r-- | src/interpolation.jl | 108 | ||||
| -rw-r--r-- | src/main.jl | 443 | ||||
| -rw-r--r-- | src/potentials.jl | 84 | ||||
| -rw-r--r-- | src/print.jl | 52 | ||||
| -rw-r--r-- | src/simpleq-Kv.jl | 119 | ||||
| -rw-r--r-- | src/simpleq-hardcore.jl | 573 | ||||
| -rw-r--r-- | src/simpleq-iteration.jl | 94 | ||||
| -rw-r--r-- | src/tools.jl | 49 |
12 files changed, 3219 insertions, 0 deletions
diff --git a/src/anyeq.jl b/src/anyeq.jl new file mode 100644 index 0000000..8459cb0 --- /dev/null +++ b/src/anyeq.jl @@ -0,0 +1,949 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# interpolation +@everywhere mutable struct Anyeq_approx + aK::Float64 + bK::Float64 + gK::Float64 + aL1::Float64 + bL1::Float64 + aL2::Float64 + bL2::Float64 + gL2::Float64 + aL3::Float64 + bL3::Float64 + gL3::Float64 +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) + # 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 + + # 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] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + + @printf("% .15e % .15e\n",E,error) +end + +# compute energy as a function of rho +function anyeq_energy_rho(rhos,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) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # 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)) + + ## 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 + + 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 + + for j in 1:length(rhos) + @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[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) + # 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 + + # compute initial guess from medeq + u0s=anyeq_init_medeq([rhos[1]],N,J,k,a0,v,maxiter,tolerance) + u=u0s[1] + + for j in 1:length(rhos) + progress(j,length(rhos),10000) + # run the task + # init Newton from previous rho + (u,E,error)=anyeq_hatu(u,P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + @printf("% .15e % .15e % .15e\n",rhos[j],E,error) + # abort when the error gets too big + if error>tolerance + break + end + end +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) + # 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 + + # compute initial guess from medeq + u0s=anyeq_init_medeq([rhos[length(rhos)]],N,J,k,a0,v,maxiter,tolerance) + u=u0s[1] + + for j in 1:length(rhos) + progress(j,length(rhos),10000) + # run the task + # init Newton from previous rho + (u,E,error)=anyeq_hatu(u,P,N,J,rhos[length(rhos)+1-j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,approx) + @printf("% .15e % .15e % .15e\n",rhos[length(rhos)+1-j],real(E),error) + # abort when the error gets too big + if error>tolerance + break + end + end +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) + # 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 + # 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] + + (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 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],u[(J-1-zeta)*N+j]) + end + end +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) + # 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 + + # 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] + + (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)) + 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) + # 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 + + # 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] + + (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 eta + eta=anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + @printf("% .15e % .15e\n",eta,error) +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) + end + for j in 1:length(rhos) + append!(work[(j-1)%nw+1],j) + end + + # 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 + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,maxiter,tolerance) + + # compute u + us=Array{Array{Float64,1}}(undef,length(rhos)) + errs=Array{Float64,1}(undef,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 + (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) + end + end + + # compute eta + etas=Array{Float64}(undef,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 + etas[j]=anyeq_eta(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) + @printf("% .15e % .15e % .15e\n",rhos[j],etas[j],errs[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) + else + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + end + + # 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] + + (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 + M=anyeq_momentum(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,approx) + + 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]) + 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) + # 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 + + # 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] + + # 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) + end +end + +# compute Abar +function anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + if approx.bL3==0. + return [] + 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 + + 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,N*J,10000) + end + # run the task + out[j]=remotecall_fetch(barAmat,workers()[p],k[j],weights,taus,T,P,N,J,2,2,2,2,2) + end + end + + return out +end + +# initialize computation +@everywhere function anyeq_init(taus,P,N,J,v) + # Gauss-Legendre weights + weights=gausslegendre(N) + + # initialize vectors V,k + V=Array{Float64}(undef,J*N) + k=Array{Float64}(undef,J*N) + for zeta in 0:J-1 + for j in 1:N + xj=weights[1][j] + # set kj + k[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1]) + # set v + V[zeta*N+j]=v(k[zeta*N+j]) + end + end + # potential at 0 + V0=v(0) + + # initialize matrix A + T=chebyshev_polynomials(P) + A=Amat(k,weights,taus,T,P,N,J,2,2) + + # compute Upsilon + # Upsilonmat does not use splines, so increase precision + weights_plus=gausslegendre(N*J) + Upsilon=Upsilonmat(k,v,weights_plus) + Upsilon0=Upsilon0mat(k,v,weights_plus) + + return(weights,T,k,V,V0,A,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)) + + weights_medeq=gausslegendre(N*J) + + (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 + print(stderr,"warning: computation of initial Ansatz failed for rho=",rhos[j],"\n") + end + end + + return u0s +end + +# interpolate the solution of medeq to an input for anyeq +@everywhere function easyeq_to_anyeq(u_simple,weights,k,N,J) + # 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)) + for j in 1:length(u_simple) + xj=weights[1][j] + k_s[length(u_simple)-j+1]=(1-xj)/(1+xj) + u_s[length(u_simple)-j+1]=u_simple[j] + end + + # initialize U + u=zeros(Float64,J*N) + for zeta in 0:J-1 + for j in 1:N + u[zeta*N+j]=linear_interpolation(k[zeta*N+j],k_s,u_s) + end + end + + return u +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) + # init + # rescale by rho (that's how u is defined) + u=rho*u0 + + # quantify relative error + error=-1. + + # 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) + + error=norm(new-u)/norm(u) + if(error<tolerance) + u=new + break + end + + u=new + end + + energy=rho/2*V0-avg_v_chebyshev(u,Upsilon0,k,taus,weights,N,J)/2 + return(u/rho,energy,error) +end + + +# save Abar +function anyeq_save_Abar(taus,P,N,J,v,approx) + # initialize vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # init Abar + Abar=anyeq_Abar_multithread(k,weights,taus,T,P,N,J,approx) + + # print params + @printf("## P=%d N=%d J=%d\n",P,N,J) + + for i in 1:N*J + for j1 in 1:(P+1)*J + for j2 in 1:(P+1)*J + for j3 in 1:(P+1)*J + for j4 in 1:(P+1)*J + for j5 in 1:(P+1)*J + @printf("% .15e\n",Abar[i][j1,j2,j3,j4,j5]) + end + end + end + end + end + end +end + +# read Abar +function anyeq_read_Abar(savefile,P,N,J) + # open file + ff=open(savefile,"r") + # read all lines + lines=readlines(ff) + close(ff) + + # init + Abar=Array{Array{Float64,5}}(undef,N*J) + for i in 1:N*J + Abar[i]=Array{Float64,5}(undef,(P+1)*J,(P+1)*J,(P+1)*J,(P+1)*J,(P+1)*J) + end + + i=1 + j1=1 + j2=1 + j3=1 + j4=1 + j5=0 + for l in 1:length(lines) + # drop comments + if lines[l]!="" && lines[l][1]!='#' + # increment counters + if j5<(P+1)*J + j5+=1 + else + j5=1 + if j4<(P+1)*J + j4+=1 + else + j4=1 + if j3<(P+1)*J + j3+=1 + else + j3=1 + if j2<(P+1)*J + j2+=1 + else + j2=1 + if j1<(P+1)*J + j1+=1 + else + j1=1 + if i<N*J + i+=1 + else + print(stderr,"error: too many lines in savefile\n") + exit() + end + end + end + end + end + end + + Abar[i][j1,j2,j3,j4,j5]=parse(Float64,lines[l]) + + end + end + + return Abar + +end + + + +# Xi +# takes the vector of kj's and xn's as input +@everywhere function anyeq_Xi(U,X,Y) + 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) + out=Array{Float64,2}(undef,N*J,N*J) + for zetapp in 0:J-1 + for n in 1:N + one=zeros(Int64,J*N) + one[zetapp*N+n]=1 + + # Chebyshev expansion of U + FU=chebyshev(U,taus,weights,P,N,J,2) + + dS=-conv_one_v_chebyshev(n,zetapp,Upsilon,k,taus,weights,N,J)/rho + + dE=-avg_one_v_chebyshev(n,zetapp,Upsilon0,k,taus,weights,N)/rho + + UU=conv_chebyshev(FU,FU,A) + + dII=zeros(Float64,N*J) + if approx.gL2!=0. + if approx.bL2!=0. + dII+=approx.gL2*approx.bL2*(conv_one_chebyshev(n,zetapp,chebyshev(U.*S,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(one.*S+U.*dS,taus,weights,P,N,J,2),A)/rho) + end + if approx.bL2!=1. + dII+=approx.gL2*(1-approx.bL2)*(dE/rho*UU+2*E/rho*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2)) + end + end + + dJJ=zeros(Float64,J*N) + if approx.gL3!=0. + if approx.bL3!=0. + FS=chebyshev(S,taus,weights,P,N,J,2) + dFU=chebyshev(one,taus,weights,P,N,J,2) + dFS=chebyshev(dS,taus,weights,P,N,J,2) + dJJ+=approx.gL3*approx.bL3*(4*double_conv_S_chebyshev(FU,FU,FU,dFU,FS,Abar)+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+4*E*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2).*UU/rho^2) + end + end + + dG=zeros(Float64,N*J) + if approx.aK!=0. && approx.gK!=0. + if approx.bK!=0. + dG+=approx.aK*approx.gK*approx.bK*(conv_one_chebyshev(n,zetapp,chebyshev(2*S.*U,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*S.*one+2*dS.*U,taus,weights,P,N,J,2),A)/rho) + end + if approx.bK!=1. + dG+=approx.aK*approx.gK*(1-approx.bK)*(2*dE*UU/rho+4*E*conv_one_chebyshev(n,zetapp,FU,A,taus,weights,P,N,J,2)/rho) + end + end + if approx.aL1!=0. + if approx.bL1!=0. + dG-=approx.aL1*approx.bL1*(conv_one_chebyshev(n,zetapp,chebyshev(S.*(U.^2),taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*S.*U.*one+dS.*(U.^2),taus,weights,P,N,J,2),A)/rho) + end + if approx.bL1!=1. + dG-=approx.aL1*(1-approx.bL1)*(E/rho*conv_one_chebyshev(n,zetapp,chebyshev((U.^2),taus,weights,P,N,J,2),taus,A,weights,P,N,J,2)+conv_chebyshev(FU,chebyshev(2*E*U.*one+dE*(U.^2),taus,weights,P,N,J,2),A)/rho) + end + end + if approx.aL2!=0. && approx.gL2!=0. + dG+=approx.aL2*(conv_one_chebyshev(n,zetapp,chebyshev(2*II.*U,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(2*dII.*U+2*II.*one,taus,weights,P,N,J,2),A)/rho) + end + if approx.aL3!=0. && approx.gL3!=0. + dG-=approx.aL3*(conv_one_chebyshev(n,zetapp,chebyshev(JJ/2,taus,weights,P,N,J,2),A,taus,weights,P,N,J,2)/rho+conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho) + end + + dsK=zeros(Float64,N*J) + if approx.gK!=0. + if approx.bK!=0. + dsK+=approx.gK*approx.bK*dS + end + if approx.bK!=1. + dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J) + end + else + end + dsL1=zeros(Float64,N*J) + if approx.bL1!=0. + dsL1+=approx.bL1*dS + end + if approx.bL1!=1. + dsL1+=(1-approx.bL1)*dE*ones(Float64,N*J) + end + + dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1 + dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1 + + out[:,zetapp*N+n]=one+((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) + end + end + return out +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) + # Chebyshev expansion of U + FU=chebyshev(U,taus,weights,P,N,J,2) + + S=v-conv_v_chebyshev(U,Upsilon,k,taus,weights,N,J)/rho + E=v0-avg_v_chebyshev(U,Upsilon0,k,taus,weights,N,J)/rho + + UU=conv_chebyshev(FU,FU,A) + + II=zeros(Float64,N*J) + if approx.gL2!=0. + if approx.bL2!=0. + II+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*S,taus,weights,P,N,J,2),A)/rho + end + if approx.bL2!=1. + II+=approx.gL2*(1-approx,bL2)*E/rho*UU + end + end + + JJ=zeros(Float64,N*J) + if approx.gL3!=0. + if approx.bL3!=0. + FS=chebyshev(S,taus,weights,P,N,J,2) + JJ+=approx.gL3*approx.bL3*double_conv_S_chebyshev(FU,FU,FU,FU,FS,Abar) + end + if approx.bL3!=1. + JJ+=approx.gL3*(1-approx.bL3)*E*(UU/rho).^2 + end + end + + G=zeros(Float64,N*J) + if approx.aK!=0. && approx.gK!=0. + if approx.bK!=0. + G+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*S.*U,taus,weights,P,N,J,2),A)/rho + end + if approx.bK!=1. + G+=approx.aK*approx.gK*(1-approx.bK)*2*E*UU/rho + end + end + if approx.aL1!=0. + if approx.bL1!=0. + G-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(S.*(U.^2),taus,weights,P,N,J,2),A)/rho + end + if approx.bL1!=1. + G-=approx.aL1*(1-approx.bL1)*E/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A) + end + end + if approx.aL2!=0. && approx.gL2!=0. + G+=approx.aL2*conv_chebyshev(FU,chebyshev(2*II.*U,taus,weights,P,N,J,2),A)/rho + end + if approx.aL3!=0 && approx.gL3!=0. + G-=approx.aL3*conv_chebyshev(FU,chebyshev(JJ/2,taus,weights,P,N,J,2),A)/rho + end + + sK=zeros(Float64,N*J) + if approx.gK!=0. + if approx.bK!=0. + sK+=approx.gK*approx.bK*S + end + if approx.bK!=1. + sK+=approx.gK*(1-approx.bK)*E*ones(Float64,N*J) + end + end + + sL1=zeros(Float64,N*J) + if approx.bL1!=0. + sL1+=approx.bL1*S + end + if approx.bL1!=1. + sL1+=(1-approx.bL1)*E*ones(Float64,N*J) + end + + X=(k.^2/2+rho*(sK-sL1+II))./sL1/rho + Y=(S-sL1+G+JJ/2)./sL1 + + return(S,E,II,JJ,X,Y,sL1,sK,G) +end + +# condensate fraction +@everywhere function anyeq_eta(u,P,N,J,rho,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,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 + 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] + + # compute M + M[i]=-avg_v_chebyshev(du,Upsilon0,k,taus,weights,N,J)/2 + end + + return M +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) + # initialize dV + dV=Array{Float64}(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 + end + dV0=1. + + # 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) + + 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) + # 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 +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) + # Chebyshev expansion of U + FU=chebyshev(U,taus,weights,P,N,J,2) + + dS=dv-conv_v_chebyshev(U,dUpsilon,k,taus,weights,N,J)/rho + dE=dv0-avg_v_chebyshev(U,dUpsilon0,k,taus,weights,N,J)/rho + + UU=conv_chebyshev(FU,FU,A) + + dII=zeros(Float64,N*J) + if approx.gL2!=0. + if approx.bL2!=0. + dII+=approx.gL2*approx.bL2*conv_chebyshev(FU,chebyshev(U.*dS,taus,weights,P,N,J,2),A)/rho + end + if approx.bL2!=1. + dII+=approx.gL2*(1-approx,bL2)*dE/rho*UU + end + end + + dJJ=zeros(Float64,J*N) + if approx.gL3!=0. + if approx.bL3!=0. + dFS=chebyshev(dS,taus,weights,P,N,J,2) + 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 + end + end + + dG=zeros(Float64,N*J) + if approx.aK!=0. && approx.gK!=0. + if approx.bK!=0. + dG+=approx.aK*approx.gK*approx.bK*conv_chebyshev(FU,chebyshev(2*dS.*U,taus,weights,P,N,J,2),A)/rho + end + if approx.bK!=1. + dG+=approx.aK*approx.gK*(1-approx.bK)*2*dE*UU/rho + end + end + if approx.aL1!=0. + if approx.bL1!=0. + dG-=approx.aL1*approx.bL1*conv_chebyshev(FU,chebyshev(dS.*(U.^2),taus,weights,P,N,J,2),A)/rho + end + if approx.bL1!=1. + dG-=approx.aL1*(1-approx.bL1)*dE/rho*conv_chebyshev(FU,chebyshev((U.^2),taus,weights,P,N,J,2),A) + end + end + if approx.aL2!=0. && approx.gL2!=0. + dG+=approx.aL2*conv_chebyshev(FU,chebyshev(2*dII.*U,taus,weights,P,N,J,2),A)/rho + end + if approx.aL3!=0. && approx.gL3!=0. + dG-=approx.aL3*conv_chebyshev(FU,chebyshev(dJJ/2,taus,weights,P,N,J,2),A)/rho + end + + dsK=zeros(Float64,N*J) + if approx.gK!=0. + if approx.bK!=0. + dsK+=approx.gK*approx.bK*dS + end + if approx.bK!=1. + dsK+=approx.gK*(1-approx.bK)*dE*ones(N*J) + end + end + dsL1=zeros(Float64,N*J) + if approx.bL1!=0. + dsL1+=approx.bL1*dS + end + if approx.bL1!=1. + dsL1+=(1-approx.bL1)*dE*ones(N*J) + end + + dX=(dsK-dsL1+dII)./sL1-X./sL1.*dsL1 + dY=(dS-dsL1+dG+dJJ/2)./sL1-Y./sL1.*dsL1 + + 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 diff --git a/src/chebyshev.jl b/src/chebyshev.jl new file mode 100644 index 0000000..28c8f1f --- /dev/null +++ b/src/chebyshev.jl @@ -0,0 +1,279 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# Chebyshev expansion +@everywhere function chebyshev(a,taus,weights,P,N,J,nu) + out=zeros(Float64,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for j in 1:N + out[zeta*(P+1)+n+1]+=(2-(n==0 ? 1 : 0))/2*weights[2][j]*cos(n*pi*(1+weights[1][j])/2)*a[zeta*N+j]/(1-(taus[zeta+2]-taus[zeta+1])/2*sin(pi*weights[1][j]/2)+(taus[zeta+2]+taus[zeta+1])/2)^nu + end + end + end + + return out +end + +# evaluate function from Chebyshev expansion +@everywhere function chebyshev_eval(Fa,x,taus,chebyshev,P,J,nu) + # change variable + tau=(1-x)/(1+x) + + out=0. + for zeta in 0:J-1 + # check that the spline is right + if tau<taus[zeta+2] && tau>=taus[zeta+1] + for n in 0:P + out+=Fa[zeta*(P+1)+n+1]*chebyshev[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])) + end + break + end + end + + return (1+tau)^nu*out +end + + +# convolution +# input the Chebyshev expansion of a and b, as well as the A matrix +@everywhere function conv_chebyshev(Fa,Fb,A) + out=zeros(Float64,length(A)) + for i in 1:length(A) + out[i]=dot(Fa,A[i]*Fb) + end + return out +end + +# <ab> +@everywhere function avg_chebyshev(Fa,Fb,A0) + return dot(Fa,A0*Fb) +end + +# 1_n * a +@everywhere function conv_one_chebyshev(n,zetapp,Fa,A,taus,weights,P,N,J,nu1) + out=zeros(Float64,N*J) + for m in 1:N*J + for l in 0:P + for p in 1:J*(P+1) + out[m]+=(2-(l==0 ? 1 : 0))/2*weights[2][n]*cos(l*pi*(1+weights[1][n])/2)/(1-(taus[zetapp+2]-taus[zetapp+1])/2*sin(pi*weights[1][n]/2)+(taus[zetapp+2]+taus[zetapp+1])/2)^nu1*A[m][zetapp*(P+1)+l+1,p]*Fa[p] + end + end + end + return out +end +# a * v +@everywhere function conv_v_chebyshev(a,Upsilon,k,taus,weights,N,J) + out=zeros(Float64,J*N) + for i in 1:J*N + for zetap in 0:J-1 + for j in 1:N + xj=weights[1][j] + out[i]+=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][j]*cos(pi*xj/2)*(1+k[zetap*N+j])^2*k[zetap*N+j]*a[zetap*N+j]*Upsilon[zetap*N+j][i] + end + end + end + return out +end +@everywhere function conv_one_v_chebyshev(n,zetap,Upsilon,k,taus,weights,N,J) + out=zeros(Float64,J*N) + xj=weights[1][n] + for i in 1:J*N + out[i]=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon[zetap*N+n][i] + end + return out +end + +# <av> +@everywhere function avg_v_chebyshev(a,Upsilon0,k,taus,weights,N,J) + out=0. + for zetap in 0:J-1 + for j in 1:N + xj=weights[1][j] + out+=(taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][j]*cos(pi*xj/2)*(1+k[zetap*N+j])^2*k[zetap*N+j]*a[zetap*N+j]*Upsilon0[zetap*N+j] + end + end + return out +end +# <1_nv> +@everywhere function avg_one_v_chebyshev(n,zetap,Upsilon0,k,taus,weights,N) + xj=weights[1][n] + return (taus[zetap+2]-taus[zetap+1])/(32*pi)*weights[2][n]*cos(pi*xj/2)*(1+k[zetap*N+n])^2*k[zetap*N+n]*Upsilon0[zetap*N+n] +end + +# compute \int dq dxi u1(k-xi)u2(q)u3(xi)u4(k-q)u5(xi-q) +@everywhere function double_conv_S_chebyshev(FU1,FU2,FU3,FU4,FU5,Abar) + out=zeros(Float64,length(Abar)) + for i in 1:length(Abar) + for j1 in 1:length(FU1) + for j2 in 1:length(FU2) + for j3 in 1:length(FU3) + for j4 in 1:length(FU4) + for j5 in 1:length(FU5) + out[i]+=Abar[i][j1,j2,j3,j4,j5]*FU1[j1]*FU2[j2]*FU3[j3]*FU4[j4]*FU5[j5] + end + end + end + end + end + end + return out +end + + +# compute A +@everywhere function Amat(k,weights,taus,T,P,N,J,nua,nub) + out=Array{Array{Float64,2},1}(undef,J*N) + for i in 1:J*N + out[i]=zeros(Float64,J*(P+1),J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=1/(pi^2*k[i])*integrate_legendre(tau->(1-tau)/(1+tau)^(3-nua)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))*(alpham(k[i],tau)>taus[zetap+2] || alphap(k[i],tau)<taus[zetap+1] ? 0. : integrate_legendre(sigma->(1-sigma)/(1+sigma)^(3-nub)*T[m+1]((2*sigma-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1],alpham(k[i],tau)),min(taus[zetap+2],alphap(k[i],tau)),weights)),taus[zeta+1],taus[zeta+2],weights) + end + end + end + end + end + + return out +end + +# compute Upsilon +@everywhere function Upsilonmat(k,v,weights) + out=Array{Array{Float64,1},1}(undef,length(k)) + for i in 1:length(k) + out[i]=Array{Float64,1}(undef,length(k)) + for j in 1:length(k) + out[i][j]=integrate_legendre(s->s*v(s)/k[j],abs(k[j]-k[i]),k[j]+k[i],weights) + end + end + return out +end +@everywhere function Upsilon0mat(k,v,weights) + out=Array{Float64,1}(undef,length(k)) + for j in 1:length(k) + out[j]=2*k[j]*v(k[j]) + end + return out +end + +# alpha_- +@everywhere function alpham(k,t) + return (1-k-(1-t)/(1+t))/(1+k+(1-t)/(1+t)) +end +# alpha_+ +@everywhere function alphap(k,t) + return (1-abs(k-(1-t)/(1+t)))/(1+abs(k-(1-t)/(1+t))) +end + + +# compute \bar A +@everywhere function barAmat(k,weights,taus,T,P,N,J,nu1,nu2,nu3,nu4,nu5) + out=zeros(Float64,J*(P+1),J*(P+1),J*(P+1),J*(P+1),J*(P+1)) + for zeta1 in 0:J-1 + for n1 in 0:P + for zeta2 in 0:J-1 + for n2 in 0:P + for zeta3 in 0:J-1 + for n3 in 0:P + for zeta4 in 0:J-1 + for n4 in 0:P + for zeta5 in 0:J-1 + for n5 in 0:P + out[zeta1*(P+1)+n1+1, + zeta2*(P+1)+n2+1, + zeta3*(P+1)+n3+1, + zeta4*(P+1)+n4+1, + zeta5*(P+1)+n5+1]=1/((2*pi)^5*k^2)*integrate_legendre(tau->barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5),taus[zeta1+1],taus[zeta1+2],weights) + end + end + end + end + end + end + end + end + end + end + + return out +end +@everywhere function barAmat_int1(tau,k,taus,T,weights,nu1,nu2,nu3,nu4,nu5,zeta1,zeta2,zeta3,zeta4,zeta5,n1,n2,n3,n4,n5) + if(alpham(k,tau)<taus[zeta2+2] && alphap(k,tau)>taus[zeta2+1]) + return 2*(1-tau)/(1+tau)^(3-nu1)*T[n1+1]((2*tau-(taus[zeta1+1]+taus[zeta1+2]))/(taus[zeta1+2]-taus[zeta1+1]))*integrate_legendre(sigma->barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5),max(taus[zeta2+1],alpham(k,tau)),min(taus[zeta2+2],alphap(k,tau)),weights) + else + return 0. + end +end +@everywhere function barAmat_int2(tau,sigma,k,taus,T,weights,nu2,nu3,nu4,nu5,zeta2,zeta3,zeta4,zeta5,n2,n3,n4,n5) + return 2*(1-sigma)/(1+sigma)^(3-nu2)*T[n2+1]((2*sigma-(taus[zeta2+1]+taus[zeta2+2]))/(taus[zeta2+2]-taus[zeta2+1]))*integrate_legendre(taup->barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5),taus[zeta3+1],taus[zeta3+2],weights) +end +@everywhere function barAmat_int3(tau,sigma,taup,k,taus,T,weights,nu3,nu4,nu5,zeta3,zeta4,zeta5,n3,n4,n5) + if(alpham(k,taup)<taus[zeta4+2] && alphap(k,taup)>taus[zeta4+1]) + return 2*(1-taup)/(1+taup)^(3-nu3)*T[n3+1]((2*taup-(taus[zeta3+1]+taus[zeta3+2]))/(taus[zeta3+2]-taus[zeta3+1]))*integrate_legendre(sigmap->barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5),max(taus[zeta4+1],alpham(k,taup)),min(taus[zeta4+2],alphap(k,taup)),weights) + else + return 0. + end +end +@everywhere function barAmat_int4(tau,sigma,taup,sigmap,k,taus,T,weights,nu4,nu5,zeta4,zeta5,n4,n5) + return 2*(1-sigmap)/(1+sigmap)^(3-nu4)*T[n4+1]((2*sigma-(taus[zeta4+1]+taus[zeta4+2]))/(taus[zeta4+2]-taus[zeta4+1]))*integrate_legendre(theta->barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5),0.,2*pi,weights) +end +@everywhere function barAmat_int5(tau,sigma,taup,sigmap,theta,k,taus,T,weights,nu5,zeta5,n5) + R=barAmat_R((1-sigma)/(1+sigma),(1-tau)/(1+tau),(1-sigmap)/(1+sigmap),(1-taup)/(1+taup),theta,k) + if((1-R)/(1+R)<taus[zeta5+2] && (1-R)/(1+R)>taus[zeta5+1]) + return (2/(2+R))^nu5*T[n5+1]((2*(1-R)/(1+R)-(taus[zeta5+1]+taus[zeta5+2]))/(taus[zeta5+2]-taus[zeta5+1])) + else + return 0. + end +end +# R(s,t,s',t,theta,k) +@everywhere function barAmat_R(s,t,sp,tp,theta,k) + return sqrt(k^2*(s^2+t^2+sp^2+tp^2)-k^4-(s^2-t^2)*(sp^2-tp^2)-sqrt((4*k^2*s^2-(k^2+s^2-t^2)^2)*(4*k^2*sp^2-(k^2+sp^2-tp^2)^2))*cos(theta))/(sqrt(2)*k) +end + +# compute Chebyshev polynomials +@everywhere function chebyshev_polynomials(P) + T=Array{Polynomial}(undef,P+1) + T[1]=Polynomial([1]) + T[2]=Polynomial([0,1]) + for n in 1:P-1 + # T_n + T[n+2]=2*T[2]*T[n+1]-T[n] + end + + return T +end + +# compute \int f*u dk/(2*pi)^3 +@everywhere function integrate_f_chebyshev(f,u,k,taus,weights,N,J) + out=0. + for zeta in 0:J-1 + for i in 1:N + out+=(taus[zeta+2]-taus[zeta+1])/(16*pi)*weights[2][i]*cos(pi*weights[1][i]/2)*(1+k[zeta*N+i])^2*k[zeta*N+i]^2*u[zeta*N+i]*f(k[zeta*N+i]) + end + end + return out +end + +@everywhere function inverse_fourier_chebyshev(u,x,k,taus,weights,N,J) + out=0. + for zeta in 0:J-1 + for j in 1:N + out+=(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 out +end diff --git a/src/easyeq.jl b/src/easyeq.jl new file mode 100644 index 0000000..0bde3ab --- /dev/null +++ b/src/easyeq.jl @@ -0,0 +1,411 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# interpolation +@everywhere mutable struct Easyeq_approx + bK::Float64 + bL::Float64 +end + +# compute energy +function easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute initial guess from previous rho + (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) + for j in 2:nlrho_init + rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) + (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) + end + + # print energy + @printf("% .15e % .15e\n",real(E),err) +end + +# compute energy as a function of rho +function easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,approx) + # compute gaussian quadrature weights + weights=gausslegendre(order) + # init u + u=easyeq_init_u(a0,order,weights) + + for j in 1:length(rhos) + # compute u (init newton with previously computed u) + (u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx) + + @printf("% .15e % .15e % .15e\n",rhos[j],real(E),err) + + end +end + +# compute u(k) +function easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx) + weights=gausslegendre(order) + + # compute initial guess from previous rho + (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) + for j in 2:nlrho_init + rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) + (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) + end + + for i in 1:order + k=(1-weights[1][i])/(1+weights[1][i]) + @printf("% .15e % .15e\n",k,real(u[i])) + end +end + +# compute u(x) +function easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx) + weights=gausslegendre(order) + + # compute initial guess from previous rho + (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) + for j in 2:nlrho_init + rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) + (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) + end + + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + @printf("% .15e % .15e\n",x,real(easyeq_u_x(x,u,weights))) + end +end + +# compute 2u(x)-rho u*u(x) +function easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,approx) + weights=gausslegendre(order) + + # compute initial guess from previous rho + (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) + for j in 2:nlrho_init + rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) + (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) + end + + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + @printf("% .15e % .15e\n",x,real(easyeq_u_x(x,2*u-rho*u.*u,weights))) + end +end + +# condensate fraction +function easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,approx) + # compute gaussian quadrature weights + weights=gausslegendre(order) + + # compute initial guess from previous rho + (u,E,err)=easyeq_hatu(easyeq_init_u(a0,order,weights),order,(10.)^minlrho_init,v,maxiter,tolerance,weights,approx) + for j in 2:nlrho_init + rho_tmp=10^(minlrho_init+(log10(rho)-minlrho_init)*(j-1)/(nlrho_init-1)) + (u,E,err)=easyeq_hatu(u,order,rho_tmp,v,maxiter,tolerance,weights,approx) + end + + # compute eta + eta=easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx) + + # print energy + @printf("% .15e % .15e\n",eta,err) +end + +# condensate fraction as a function of rho +function easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,approx) + weights=gausslegendre(order) + # init u + u=easyeq_init_u(a0,order,weights) + + for j in 1:length(rhos) + # compute u (init newton with previously computed u) + (u,E,err)=easyeq_hatu(u,order,rhos[j],v,maxiter,tolerance,weights,approx) + + # compute eta + eta=easyeq_eta(u,order,rhos[j],v,maxiter,tolerance,weights,approx) + + @printf("% .15e % .15e % .15e\n",rhos[j],eta,err) + end +end + + +# initialize u +@everywhere function easyeq_init_u(a0,order,weights) + u=zeros(Float64,order) + for j in 1:order + # transformed k + k=(1-weights[1][j])/(1+weights[1][j]) + u[j]=4*pi*a0/k^2 + end + + return u +end + +# \hat u(k) computed using Newton +@everywhere function easyeq_hatu(u0,order,rho,v,maxiter,tolerance,weights,approx) + # initialize V and Eta + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + # init u + u=rho*u0 + + # iterate + err=Inf + for i in 1:maxiter-1 + new=u-inv(easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx) + + err=norm(new-u)/norm(u) + if(err<tolerance) + u=new + break + end + u=new + end + + return (u/rho,easyeq_en(u,V0,Eta0,rho,weights)*rho/2,err) +end + +# \Eta +@everywhere function easyeq_H(x,t,weights,v) + return (x>t ? 2*t/x : 2)* integrate_legendre(y->2*pi*((x+t)*y+abs(x-t)*(1-y))*v((x+t)*y+abs(x-t)*(1-y)),0,1,weights) +end + +# initialize V +@everywhere function easyeq_init_v(weights,v) + order=length(weights[1]) + V=Array{Float64}(undef,order) + V0=v(0) + for i in 1:order + k=(1-weights[1][i])/(1+weights[1][i]) + V[i]=v(k) + end + return(V,V0) +end + +# initialize Eta +@everywhere function easyeq_init_H(weights,v) + order=length(weights[1]) + Eta=Array{Array{Float64}}(undef,order) + Eta0=Array{Float64}(undef,order) + for i in 1:order + k=(1-weights[1][i])/(1+weights[1][i]) + Eta[i]=Array{Float64}(undef,order) + for j in 1:order + y=(weights[1][j]+1)/2 + Eta[i][j]=easyeq_H(k,(1-y)/y,weights,v) + end + y=(weights[1][i]+1)/2 + Eta0[i]=easyeq_H(0,(1-y)/y,weights,v) + end + return(Eta,Eta0) +end + +# Xi(u) +@everywhere function easyeq_Xi(u,V,V0,Eta,Eta0,weights,rho,approx) + order=length(weights[1]) + + # init + out=zeros(Float64,order) + + # compute E before running the loop + E=easyeq_en(u,V0,Eta0,rho,weights) + + for i in 1:order + # k_i + k=(1-weights[1][i])/(1+weights[1][i]) + # S_i + S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights) + + # A_K,i + A=0. + if approx.bK!=0. + A+=approx.bK*S + end + if approx.bK!=1. + A+=(1-approx.bK)*E + end + + # T + if approx.bK==1. + T=1. + else + T=S/A + end + + # B + if approx.bK==approx.bL + B=1. + else + B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E)) + end + + # X_i + X=k^2/(2*A*rho) + + # U_i + out[i]=u[i]-T/(2*(X+1))*Phi(B*T/(X+1)^2) + end + + return out +end + +# derivative of Xi +@everywhere function easyeq_dXi(u,V,V0,Eta,Eta0,weights,rho,approx) + order=length(weights[1]) + + # init + out=zeros(Float64,order,order) + + # compute E before the loop + E=easyeq_en(u,V0,Eta0,rho,weights) + + for i in 1:order + # k_i + k=(1-weights[1][i])/(1+weights[1][i]) + # S_i + S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights) + + # A_K,i + A=0. + if approx.bK!=0. + A+=approx.bK*S + end + if approx.bK!=1. + A+=(1-approx.bK)*E + end + + # T + if approx.bK==1. + T=1. + else + T=S/A + end + + # B + if approx.bK==approx.bL + B=1. + else + B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E)) + end + + # X_i + X=k^2/(2*A*rho) + + for j in 1:order + y=(weights[1][j]+1)/2 + dS=-1/rho*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3)*weights[2][j] + dE=-1/rho*(1-y)*Eta0[j]/(2*(2*pi)^3*y^3)*weights[2][j] + + # dA + dA=0. + if approx.bK!=0. + dA+=approx.bK*dS + end + if approx.bK!=1. + dA+=(1-approx.bK)*dE + end + + # dT + if approx.bK==1. + dT=0. + else + dT=(1-approx.bK)*(E*dS-S*dE)/A^2 + end + + # dB + if approx.bK==approx.bL + dB=0. + else + dB=(approx.bL*(1-approx.bK)-approx.bK*(1-approx.bL))*(E*dS-S*dE)/(approx.bK*S+(1-approx.bK*E))^2 + end + + dX=-k^2/(2*A^2*rho)*dA + + out[i,j]=(i==j ? 1 : 0)-(dT-T*dX/(X+1))/(2*(X+1))*Phi(B*T/(X+1)^2)-T/(2*(X+1)^3)*(B*dT+T*dB-2*B*T*dX/(X+1))*dPhi(B*T/(X+1)^2) + end + end + + return out +end + +# derivative of Xi with respect to mu +@everywhere function easyeq_dXidmu(u,V,V0,Eta,Eta0,weights,rho,approx) + order=length(weights[1]) + + # init + out=zeros(Float64,order) + + # compute E before running the loop + E=easyeq_en(u,V0,Eta0,rho,weights) + + for i in 1:order + # k_i + k=(1-weights[1][i])/(1+weights[1][i]) + # S_i + S=V[i]-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta[i].*u,0,1,weights) + + # A_K,i + A=0. + if approx.bK!=0. + A+=approx.bK*S + end + if approx.bK!=1. + A+=(1-approx.bK)*E + end + + # T + if approx.bK==1. + T=1. + else + T=S/A + end + + # B + if approx.bK==approx.bL + B=1. + else + B=(approx.bL*S+(1-approx.bL*E))/(approx.bK*S+(1-approx.bK*E)) + end + + # X_i + X=k^2/(2*A*rho) + + out[i]=T/(2*rho*A*(X+1)^2)*Phi(B*T/(X+1)^2)+B*T^2/(rho*A*(X+1)^4)*dPhi(B*T/(X+1)^2) + end + + return out +end + +# energy +@everywhere function easyeq_en(u,V0,Eta0,rho,weights) + return V0-1/(rho*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*u,0,1,weights) +end + +# condensate fraction +@everywhere function easyeq_eta(u,order,rho,v,maxiter,tolerance,weights,approx) + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + du=-inv(easyeq_dXi(rho*u,V,V0,Eta,Eta0,weights,rho,approx))*easyeq_dXidmu(rho*u,V,V0,Eta,Eta0,weights,rho,approx) + + eta=-1/(2*(2*pi)^3)*integrate_legendre_sampled(y->(1-y)/y^3,Eta0.*du,0,1,weights) + + return eta +end + +# inverse Fourier transform +@everywhere function easyeq_u_x(x,u,weights) + order=length(weights[1]) + out=integrate_legendre_sampled(y->(1-y)/y^3*sin(x*(1-y)/y)/x/(2*pi^2),u,0,1,weights) + return out +end diff --git a/src/integration.jl b/src/integration.jl new file mode 100644 index 0000000..9be4641 --- /dev/null +++ b/src/integration.jl @@ -0,0 +1,58 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# approximate \int_a^b f using Gauss-Legendre quadratures +@everywhere function integrate_legendre(f,a,b,weights) + out=0 + for i in 1:length(weights[1]) + out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2) + end + return out +end +# \int f*g where g is sampled at the Legendre nodes +@everywhere function integrate_legendre_sampled(f,g,a,b,weights) + out=0 + for i in 1:length(weights[1]) + out+=(b-a)/2*weights[2][i]*f((b-a)/2*weights[1][i]+(b+a)/2)*g[i] + end + return out +end + + +# approximate \int_a^b f/sqrt((b-x)(x-a)) using Gauss-Chebyshev quadratures +@everywhere function integrate_chebyshev(f,a,b,N) + out=0 + for i in 1:N + out=out+pi/N*f((b-a)/2*cos((2*i-1)/(2*N)*pi)+(b+a)/2) + end + return out +end + +# approximate \int_0^\infty dr f(r)*exp(-a*r) using Gauss-Chebyshev quadratures +@everywhere function integrate_laguerre(f,a,weights_gL) + out=0. + for i in 1:length(weights_gL[1]) + out+=1/a*f(weights_gL[1][i]/a)*weights_gL[2][i] + end + return out +end + +# Hann window +@everywhere function hann(x,L) + if abs(x)<L/2 + return cos(pi*x/L)^2 + else + return 0. + end +end diff --git a/src/interpolation.jl b/src/interpolation.jl new file mode 100644 index 0000000..fa3bcdb --- /dev/null +++ b/src/interpolation.jl @@ -0,0 +1,108 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# linear interpolation: given vectors x,y, compute a linear interpolation for y(x0) +# assume x is ordered +@everywhere function linear_interpolation(x0,x,y) + # if x0 is beyond all x's, then return the corresponding boundary value. + if x0>x[length(x)] + return y[length(y)] + elseif x0<x[1] + return y[1] + end + + # find bracketing interval + i=bracket(x0,x,1,length(x)) + + # interpolate + return y[i]+(y[i+1]-y[i])*(x0-x[i])/(x[i+1]-x[i]) +end +@everywhere function bracket(x0,x,a,b) + i=floor(Int64,(a+b)/2) + if x0<x[i] + return bracket(x0,x,a,i) + elseif x0>x[i+1] + return bracket(x0,x,i+1,b) + else + return i + end +end + + +# polynomial interpolation of a family of points +@everywhere function poly_interpolation(x,y) + # init for recursion + rec=Array{Polynomial{Float64}}(undef,length(x)) + for i in 1:length(x) + rec[i]=Polynomial([1.]) + end + + # compute \prod (x-x_i) + poly_interpolation_rec(rec,x,1,length(x)) + + # sum terms together + out=0. + for i in 1:length(y) + out+=rec[i]/rec[i](x[i])*y[i] + end + + return out +end +# recursive helper function +@everywhere function poly_interpolation_rec(out,x,a,b) + if a==b + return + end + # mid point + mid=floor(Int64,(a+b)/2) + # multiply left and right of mid + right=Polynomial([1.]) + for i in mid+1:b + right*=Polynomial([-x[i],1.]) + end + left=Polynomial([1.]) + for i in a:mid + left*=Polynomial([-x[i],1.]) + end + + # multiply into left and right + for i in a:mid + out[i]*=right + end + for i in mid+1:b + out[i]*=left + end + + # recurse + poly_interpolation_rec(out,x,a,mid) + poly_interpolation_rec(out,x,mid+1,b) + + return +end +## the following does the same, but has complexity N^2, the function above has N*log(N) +#@everywhere function poly_interpolation(x,y) +# out=Polynomial([0.]) +# for i in 1:length(x) +# prod=Polynomial([1.]) +# for j in 1:length(x) +# if j!=i +# prod*=Polynomial([-x[j]/(x[i]-x[j]),1/(x[i]-x[j])]) +# end +# end +# out+=prod*y[i] +# end +# +# return out +#end + diff --git a/src/main.jl b/src/main.jl new file mode 100644 index 0000000..382fe6b --- /dev/null +++ b/src/main.jl @@ -0,0 +1,443 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +using Distributed +@everywhere using FastGaussQuadrature +@everywhere using LinearAlgebra +@everywhere using Polynomials +@everywhere using Printf +@everywhere using SpecialFunctions +include("chebyshev.jl") +include("integration.jl") +include("interpolation.jl") +include("tools.jl") +include("potentials.jl") +include("print.jl") +include("easyeq.jl") +include("simpleq-Kv.jl") +include("simpleq-iteration.jl") +include("simpleq-hardcore.jl") +include("anyeq.jl") +function main() + ## defaults + + # values of rho,e, when needed + rho=1e-6 + e=1e-4 + + # incrementally initialize Newton algorithm + nlrho_init=1 + + # potential + v=k->v_exp(k,1.) + a0=a0_exp(1.) + # arguments of the potential + v_param_a=1. + v_param_b=1. + v_param_c=1. + v_param_e=1. + + # plot range when plotting in rho + minlrho=-6 + maxlrho=2 + nlrho=100 + rhos=Array{Float64}(undef,0) + # plot range when plotting in e + minle=-6 + maxle=2 + nle=100 + es=Array{Float64}(undef,0) + # plot range when plotting in x + xmin=0 + xmax=100 + nx=100 + + # cutoffs + tolerance=1e-11 + order=100 + maxiter=21 + + # for anyeq + # P + P=11 + # N + N=12 + # number of splines + J=10 + + # starting rho from which to incrementally initialize Newton algorithm + # default must be set after reading rho, if not set explicitly + minlrho_init=nothing + + # Hann window for Fourier transforms + windowL=1e3 + + # approximation for easyeq + # bK,bL + easyeq_simpleq_approx=Easyeq_approx(0.,0.) + easyeq_medeq_approx=Easyeq_approx(1.,1.) + easyeq_approx=easyeq_simpleq_approx + + # approximation for anyeq + # aK,bK,gK,aL1,bL1,aL2,bL2,gL2,aL3,bL3,gl3 + anyeq_simpleq_approx=Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.) + anyeq_bigeq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,0.,0.,0.) + anyeq_fulleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,0.,1.) + anyeq_medeq_approx=Anyeq_approx(0.,1.,1.,0.,1.,0.,0.,0.,0.,0.,0.) + anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.) + anyeq_approx=anyeq_bigeq_approx + + # read cli arguments + (params,potential,method,savefile,command)=read_args(ARGS) + + # read params + if params!="" + for param in split(params,";") + terms=split(param,"=") + if length(terms)!=2 + print(stderr,"error: could not read parameter '",param,"'.\n") + exit(-1) + end + lhs=terms[1] + rhs=terms[2] + if lhs=="rho" + rho=parse(Float64,rhs) + elseif lhs=="minlrho_init" + minlrho_init=parse(Float64,rhs) + elseif lhs=="nlrho_init" + nlrho_init=parse(Int64,rhs) + elseif lhs=="e" + e=parse(Float64,rhs) + elseif lhs=="tolerance" + tolerance=parse(Float64,rhs) + elseif lhs=="order" + order=parse(Int64,rhs) + elseif lhs=="maxiter" + maxiter=parse(Int64,rhs) + elseif lhs=="rhos" + rhos=parse_list(rhs) + elseif lhs=="minlrho" + minlrho=parse(Float64,rhs) + elseif lhs=="maxlrho" + maxlrho=parse(Float64,rhs) + elseif lhs=="nlrho" + nlrho=parse(Int64,rhs) + elseif lhs=="es" + es=parse_list(rhs) + elseif lhs=="minle" + minle=parse(Float64,rhs) + elseif lhs=="maxle" + maxle=parse(Float64,rhs) + elseif lhs=="nle" + nle=parse(Int64,rhs) + elseif lhs=="xmin" + xmin=parse(Float64,rhs) + elseif lhs=="xmax" + xmax=parse(Float64,rhs) + elseif lhs=="nx" + nx=parse(Int64,rhs) + elseif lhs=="P" + P=parse(Int64,rhs) + elseif lhs=="N" + N=parse(Int64,rhs) + elseif lhs=="J" + J=parse(Int64,rhs) + elseif lhs=="window_L" + windowL=parse(Float64,rhs) + elseif lhs=="aK" + anyeq_approx.aK=parse(Float64,rhs) + elseif lhs=="bK" + anyeq_approx.bK=parse(Float64,rhs) + easyeq_approx.bK=parse(Float64,rhs) + elseif lhs=="gK" + anyeq_approx.gK=parse(Float64,rhs) + elseif lhs=="aL1" + anyeq_approx.aL1=parse(Float64,rhs) + elseif lhs=="bL" + easyeq_approx.bL=parse(Float64,rhs) + elseif lhs=="bL1" + anyeq_approx.bL1=parse(Float64,rhs) + elseif lhs=="aL2" + anyeq_approx.aL2=parse(Float64,rhs) + elseif lhs=="bL2" + anyeq_approx.bL2=parse(Float64,rhs) + elseif lhs=="gL2" + anyeq_approx.gL2=parse(Float64,rhs) + elseif lhs=="aL3" + anyeq_approx.aL3=parse(Float64,rhs) + elseif lhs=="bL3" + anyeq_approx.bL3=parse(Float64,rhs) + elseif lhs=="gL3" + anyeq_approx.gL3=parse(Float64,rhs) + elseif lhs=="v_a" + v_param_a=parse(Float64,rhs) + elseif lhs=="v_b" + v_param_b=parse(Float64,rhs) + elseif lhs=="v_c" + v_param_c=parse(Float64,rhs) + elseif lhs=="v_e" + v_param_e=parse(Float64,rhs) + elseif lhs=="eq" + if rhs=="simpleq" + easyeq_approx=easyeq_simpleq_approx + anyeq_approx=anyeq_simpleq_approx + elseif rhs=="medeq" + easyeq_approx=easyeq_medeq_approx + anyeq_approx=anyeq_medeq_approx + elseif rhs=="bigeq" + anyeq_approx=anyeq_bigeq_approx + elseif rhs=="fulleq" + anyeq_approx=anyeq_fulleq_approx + elseif rhs=="compleq" + anyeq_approx=anyeq_compleq_approx + else + print(stderr,"error: unrecognized equation: ",rhs,"\n") + exit(-1) + end + else + print(stderr,"unrecognized parameter '",lhs,"'.\n") + exit(-1) + end + end + end + + ## read potential + if potential=="exp" + v=k->v_exp(k,v_param_a) + a0=a0_exp(v_param_a) + elseif potential=="expcry" + v=k->v_expcry(k,v_param_a,v_param_b) + a0=a0_expcry(v_param_a,v_param_b) + elseif potential=="npt" + v=k->v_npt(k) + a0=a0_npt() + elseif potential=="alg" + v=v_alg + a0=a0_alg + elseif potential=="algwell" + v=v_algwell + a0=a0_algwell + elseif potential=="exact" + v=k->v_exact(k,v_param_b,v_param_c,v_param_e) + a0=a0_exact(v_param_b,v_param_c,v_param_e) + elseif potential=="tent" + v=k->v_tent(k,v_param_a,v_param_b) + a0=a0_tent(v_param_a,v_param_b) + else + print(stderr,"unrecognized potential '",potential,"'.\n'") + exit(-1) + end + + ## set parameters + # rhos + if length(rhos)==0 + rhos=Array{Float64}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j)) + end + end + # es + if length(es)==0 + es=Array{Float64}(undef,nle) + for j in 0:nle-1 + es[j+1]=(nle==1 ? 10^minle : 10^(minle+(maxle-minle)/(nle-1)*j)) + end + end + + # default minlrho_init + if (minlrho_init==nothing) + minlrho_init=log10(rho) + end + + # splines + taus=Array{Float64}(undef,J+1) + for j in 0:J + taus[j+1]=-1+2*j/J + end + + ## run command + if method=="easyeq" + if command=="energy" + easyeq_energy(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) + # e(rho) + elseif command=="energy_rho" + easyeq_energy_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx) + # u(k) + elseif command=="uk" + easyeq_uk(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) + # u(x) + elseif command=="ux" + easyeq_ux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,easyeq_approx) + elseif command=="uux" + easyeq_uux(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,easyeq_approx) + # condensate fraction + elseif command=="condensate_fraction" + easyeq_condensate_fraction(minlrho_init,nlrho_init,order,rho,a0,v,maxiter,tolerance,easyeq_approx) + elseif command=="condensate_fraction_rho" + easyeq_condensate_fraction_rho(rhos,order,a0,v,maxiter,tolerance,easyeq_approx) + else + print(stderr,"unrecognized command '",command,"'.\n") + exit(-1) + end + elseif method=="simpleq-hardcore" + if command=="energy_rho" + simpleq_hardcore_energy_rho(rhos,taus,P,N,J,maxiter,tolerance) + elseif command=="ux" + simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance) + elseif command=="condensate_fraction_rho" + simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,maxiter,tolerance) + end + elseif method=="simpleq-iteration" + # u_n(x) using iteration + if command=="u" + simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx) + # rho(e) using iteration + elseif command=="rho_e" + simpleq_iteration_rho_e(es,order,v,maxiter) + else + print(stderr,"unrecognized command '",command,"'.\n") + exit(-1) + end + elseif method=="anyeq" + if command=="save_Abar" + anyeq_save_Abar(taus,P,N,J,v,anyeq_approx) + elseif command=="energy" + anyeq_energy(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + # e(rho) + elseif command=="energy_rho" + anyeq_energy_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + elseif command=="energy_rho_init_prevrho" + anyeq_energy_rho_init_prevrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + elseif command=="energy_rho_init_nextrho" + anyeq_energy_rho_init_nextrho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + # u(j) + elseif command=="uk" + anyeq_uk(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,anyeq_approx,savefile) + # u(j) + elseif command=="ux" + anyeq_ux(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile) + # condensate fraction + elseif command=="condensate_fraction" + anyeq_condensate_fraction(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + elseif command=="condensate_fraction_rho" + anyeq_condensate_fraction_rho(rhos,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + # momentum distribution + elseif command=="momentum_distribution" + anyeq_momentum_distribution(rho,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) + elseif command=="2pt" + anyeq_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile) + else + print(stderr,"unrecognized command: '",command,"'\n") + exit(-1) + end + elseif method=="simpleq-Kv" + if command=="2pt" + simpleq_Kv_2pt(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) + elseif command=="condensate_fraction_rho" + simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance) + elseif command=="Kv" + simpleq_Kv_Kv(minlrho_init,nlrho_init,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) + else + print(stderr,"unrecognized command: '",command,"'\n") + exit(-1) + end + else + print(stderr,"unrecognized method: '",method,"'\n") + exit(-1) + end + +end + +# parse a comma separated list as an array of Float64 +function parse_list(str) + elems=split(str,",") + out=Array{Float64}(undef,length(elems)) + for i in 1:length(elems) + out[i]=parse(Float64,elems[i]) + end + return out +end + +# read cli arguments +function read_args(ARGS) + # flag + flag="" + + # output strings + params="" + # default potential + potential="exp" + # default method + method="easyeq" + + savefile="" + command="" + + # loop over arguments + for arg in ARGS + # argument that starts with a dash + if arg[1]=='-' + # go through characters after dash + for char in arg[2:length(arg)] + + # set params + if char=='p' + # raise flag + flag="params" + elseif char=='U' + # raise flag + flag="potential" + elseif char=='M' + # raise flag + flag="method" + elseif char=='s' + # raise flag + flag="savefile" + else + print_usage() + exit(-1) + end + end + # arguments that do not start with a dash + else + if flag=="params" + params=arg + elseif flag=="potential" + potential=arg + elseif flag=="method" + method=arg + elseif flag=="savefile" + savefile=arg + else + command=arg + end + # reset flag + flag="" + end + end + + if command=="" + print_usage() + exit(-1) + end + + return (params,potential,method,savefile,command) +end + +# usage +function print_usage() + print(stderr,"usage: simplesolv [-p params] [-U potential] [-M method] [-s savefile] <command>\n") +end + +main() diff --git a/src/potentials.jl b/src/potentials.jl new file mode 100644 index 0000000..46cafc0 --- /dev/null +++ b/src/potentials.jl @@ -0,0 +1,84 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# exponential potential in 3 dimensions +@everywhere function v_exp(k,a) + return 8*pi/(1+k^2)^2*a +end +@everywhere function a0_exp(a) + if a>0. + return log(a)+2*MathConstants.eulergamma+2*besselk(0,2*sqrt(a))/besseli(0,2*sqrt(a)) + elseif a<0. + return log(-a)+2*MathConstants.eulergamma-pi*bessely(0,2*sqrt(-a))/besselj(0,2*sqrt(-a)) + else + return 0. + end +end + +# exp(-x)-a*exp(-b*x) in 3 dimensions +@everywhere function v_expcry(k,a,b) + return 8*pi*((1+k^2)^(-2)-a*b*(b^2+k^2)^(-2)) +end +@everywhere function a0_expcry(a,b) + return 1.21751642717932720441274114683413710125487579284827462 #ish +end + +# x^2*exp(-|x|) in 3 dimensions +@everywhere function v_npt(k) + return 96*pi*(1-k^2)/(1+k^2)^4 +end +@everywhere function a0_npt() + return 1. #ish +end + +# 1/(1+x^4/4) potential in 3 dimensions +@everywhere function v_alg(k) + if(k==0) + return 4*pi^2 + else + return 4*pi^2*exp(-k)*sin(k)/k + end +end +a0_alg=1. #ish + +# (1+a x^4)/(1+x^2)^4 potential in 3 dimensions +@everywhere function v_algwell(k) + a=4 + return pi^2/24*exp(-k)*(a*(k^2-9*k+15)+k^2+3*k+3) +end +a0_algwell=1. #ish + +# potential corresponding to the exact solution c/(1+b^2x^2)^2 +@everywhere function v_exact(k,b,c,e) + if k!=0 + return 48*pi^2*((18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))*exp(-sqrt(1-sqrt(c))*k/b)+(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))*exp(-sqrt(1+sqrt(c))*k/b)+(1-k/b)/2*exp(-k/b)-c*e/b^2*(3*(9-c)*k/b+8*c)/(8*(9-c)^2)*exp(-2*k/b))/k + else + return 48*pi^2*(-sqrt(1-sqrt(c))/b*(18+3*sqrt(c)-(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3+sqrt(c))^2*sqrt(c))-sqrt(1+sqrt(c))/b*(-18+3*sqrt(c)+(4-3*e/b^2)*c-(1-2*e/b^2)*c^1.5)/(4*(3-sqrt(c))^2*sqrt(c))-1/b-c*e/b^2*(27-19*c)/(8*(9-c)^2)) + end +end +@everywhere function a0_exact(b,c,e) + return 1. #ish +end + +# tent potential (convolution of soft sphere with itself): a*pi/12*(2*|x|/b-2)^2*(2*|x|/b+4) for |x|<b +@everywhere function v_tent(k,a,b) + if k!=0 + return (b/2)^3*a*(4*pi*(sin(k*b/2)-k*b/2*cos(k*b/2))/(k*b/2)^3)^2 + else + return (b/2)^3*a*(4*pi/3)^2 + end +end +@everywhere function a0_tent(a,b) + return b #ish +end diff --git a/src/print.jl b/src/print.jl new file mode 100644 index 0000000..bef1c4d --- /dev/null +++ b/src/print.jl @@ -0,0 +1,52 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# print progress +@everywhere function progress( + j::Int64, + tot::Int64, + freq::Int64 +) + if (j-1)%ceil(Int,tot/freq)==0 + if j>1 + @printf(stderr,"\r") + end + @printf(stderr,"%d/%d",j,tot) + end + if j==tot + @printf(stderr,"\r") + @printf(stderr,"%d/%d\n",j,tot) + end +end + +# print progress of two indices at once +@everywhere function progress_mat( + j1::Int64, + tot1::Int64, + j2::Int64, + tot2::Int64, + freq::Int64 +) + if ((j1-1)*tot2+j2-1)%ceil(Int,tot1*tot2/freq)==0 + if j1>1 || j2>1 + @printf(stderr,"\r") + end + @printf(stderr,"%2d/%2d, %2d/%2d",j1,tot1,j2,tot2) + end + if j1==tot1 && j2==tot2 + @printf(stderr,"\r") + @printf(stderr,"%2d/%2d, %2d/%2d\n",j1,tot1,j2,tot2) + end +end + diff --git a/src/simpleq-Kv.jl b/src/simpleq-Kv.jl new file mode 100644 index 0000000..8789656 --- /dev/null +++ b/src/simpleq-Kv.jl @@ -0,0 +1,119 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# Compute Kv=(-\Delta+v+4e(1-\rho u*))^{-1}v +function anyeq_Kv(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) + # init vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # 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] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)) + + # Kv in Fourier space + Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) + + # switch to real space + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + Kv=0. + for zeta in 0:J-1 + for j in 1:N + Kv+=(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]*Kvk[zeta*N+j]*sin(k[zeta*N+j]*x) + end + end + @printf("% .15e % .15e\n",x,Kv) + end +end + +# Compute the condensate fraction for simpleq using Kv +function simpleq_Kv_condensate_fraction(rhos,taus,P,N,J,a0,v,maxiter,tolerance) + # init vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # compute initial guess from medeq + u0s=anyeq_init_medeq(rhos,N,J,k,a0,v,0.,maxiter,tolerance) + + for j in 1:length(rhos) + (u,E,error)=anyeq_hatu(u0s[j],0.,P,N,J,rhos[j],a0,weights,k,taus,V,V0,A,Abar,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)) + + # Kv in Fourier space + Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) + + # denominator: (1-rho*\int (Kv)(2u-rho u*u)) + denom=1-rhos[j]*integrate_f_chebyshev(s->1.,Kvk.*(2*u-rhos[j]*u.^2),k,taus,weights,N,J) + + eta=rhos[j]*integrate_f_chebyshev(s->1.,Kvk.*u,k,taus,weights,N,J)/denom + + @printf("% .15e % .15e\n",rhos[j],eta) + end +end + +# Compute the two-point correlation function for simpleq using Kv +function simpleq_Kv_2pt(minlrho,nlrho,taus,P,N,J,rho,a0,v,maxiter,tolerance,xmin,xmax,nx) + # init vectors + (weights,T,k,V,V0,A,Upsilon,Upsilon0)=anyeq_init(taus,P,N,J,v) + + # 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] + + (u,E,error)=anyeq_hatu(u0,P,N,J,rho,a0,weights,k,taus,V,V0,A,nothing,Upsilon,Upsilon0,v,maxiter,tolerance,Anyeq_approx(0.,0.,1.,0.,0.,0.,0.,0.,0.,0.,0.)) + + # Kv in Fourier space + Kvk=simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) + + # denominator: (1-rho*\int (Kv)(2u-rho u*u)) + denom=1-rho*integrate_f_chebyshev(s->1.,Kvk.*(2*u-rho*u.^2),k,taus,weights,N,J) + + # Kv, u, u*Kv, u*u*Kv in real space + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + ux=inverse_fourier_chebyshev(u,x,k,taus,weights,N,J) + Kv=inverse_fourier_chebyshev(Kvk,x,k,taus,weights,N,J) + uKv=inverse_fourier_chebyshev(u.*Kvk,x,k,taus,weights,N,J) + uuKv=inverse_fourier_chebyshev(u.*u.*Kvk,x,k,taus,weights,N,J) + + # correlation in real space + Cx=rho^2*((1-ux)-((1-ux)*Kv-2*rho*uKv+rho^2*uuKv)/denom) + + @printf("% .15e % .15e\n",x,Cx) + end +end + +# Kv +function simpleq_Kv_Kvk(u,V,E,rho,Upsilon,k,taus,weights,N,J) + # (-Delta+v+4e(1-\rho u*)) in Fourier space + M=Array{Float64,2}(undef,N*J,N*J) + for zetapp in 0:J-1 + for n in 1:N + M[:,zetapp*N+n]=conv_one_v_chebyshev(n,zetapp,Upsilon,k,taus,weights,N,J) + end + end + for i in 1:J*N + M[i,i]+=k[i]^2+4*E*(1-rho*u[i]) + end + + return inv(M)*V +end diff --git a/src/simpleq-hardcore.jl b/src/simpleq-hardcore.jl new file mode 100644 index 0000000..ca64f78 --- /dev/null +++ b/src/simpleq-hardcore.jl @@ -0,0 +1,573 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# compute energy as a function of rho +function simpleq_hardcore_energy_rho(rhos,taus,P,N,J,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%nw+1],j) + end + + # initialize vectors + (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) + + # initial guess + u0s=Array{Array{Float64}}(undef,length(rhos)) + e0s=Array{Float64}(undef,length(rhos)) + for j in 1:length(rhos) + u0s[j]=Array{Float64}(undef,N*J) + for i in 1:N*J + u0s[j][i]=1/(1+r[i]^2)^2 + end + e0s[j]=pi*rhos[j] + end + + + # save result from each task + us=Array{Array{Float64}}(undef,length(rhos)) + es=Array{Float64}(undef,length(rhos)) + err=Array{Float64}(undef,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 + (us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + end + end + + for j in 1:length(rhos) + @printf("% .15e % .15e % .15e\n",rhos[j],es[j],err[j]) + end +end + +# compute u(x) +function simpleq_hardcore_ux(rho,taus,P,N,J,maxiter,tolerance) + # initialize vectors + (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) + + # initial guess + u0=Array{Float64}(undef,N*J) + for i in 1:N*J + u0[i]=1/(1+r[i]^2)^2 + end + e0=pi*rho + + (u,e,err)=simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + + for i in 1:N*J + @printf("% .15e % .15e\n",r[i],u[i]) + end +end + +# compute condensate fraction as a function of rho +function simpleq_hardcore_condensate_fraction_rho(rhos,taus,P,N,J,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%nw+1],j) + end + + # initialize vectors + (weights,weights_gL,r,T)=simpleq_hardcore_init(taus,P,N,J) + + # initial guess + u0s=Array{Array{Float64}}(undef,length(rhos)) + e0s=Array{Float64}(undef,length(rhos)) + for j in 1:length(rhos) + u0s[j]=Array{Float64}(undef,N*J) + for i in 1:N*J + u0s[j][i]=1/(1+r[i]^2)^2 + end + e0s[j]=pi*rhos[j] + end + + + # save result from each task + us=Array{Array{Float64}}(undef,length(rhos)) + es=Array{Float64}(undef,length(rhos)) + err=Array{Float64}(undef,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 + (us[j],es[j],err[j])=remotecall_fetch(simpleq_hardcore_hatu,workers()[p],u0s[j],e0s[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + end + end + + for j in 1:length(rhos) + du=-inv(simpleq_hardcore_DXi(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_dXidmu(us[j],es[j],rhos[j],r,taus,T,weights,weights_gL,P,N,J) + @printf("% .15e % .15e % .15e\n",rhos[j],du[N*J+1],err[j]) + end +end + + +# initialize computation +@everywhere function simpleq_hardcore_init(taus,P,N,J) + # Gauss-Legendre weights + weights=gausslegendre(N) + weights_gL=gausslaguerre(N) + + # r + r=Array{Float64}(undef,J*N) + for zeta in 0:J-1 + for j in 1:N + xj=weights[1][j] + # set kj + r[zeta*N+j]=(2+(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)-(taus[zeta+2]+taus[zeta+1]))/(2-(taus[zeta+2]-taus[zeta+1])*sin(pi*xj/2)+taus[zeta+2]+taus[zeta+1]) + end + end + + # Chebyshev polynomials + T=chebyshev_polynomials(P) + + return (weights,weights_gL,r,T) +end + +# compute u using chebyshev expansions +@everywhere function simpleq_hardcore_hatu(u0,e0,rho,r,taus,T,weights,weights_gL,P,N,J,maxiter,tolerance) + # init + vec=Array{Float64}(undef,J*N+1) + for i in 1:J*N + vec[i]=u0[i] + end + vec[J*N+1]=e0 + + # quantify relative error + error=Inf + + # run Newton algorithm + for i in 1:maxiter-1 + u=vec[1:J*N] + e=vec[J*N+1] + new=vec-inv(simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J))*simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + + error=norm(new-vec)/norm(new) + if(error<tolerance) + vec=new + break + end + + vec=new + end + + return(vec[1:J*N],vec[J*N+1],error) +end + +# Xi +@everywhere function simpleq_hardcore_Xi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + out=Array{Float64}(undef,J*N+1) + FU=chebyshev(u,taus,weights,P,N,J,4) + + #D's + d0=D0(e,rho,r,weights,N,J) + d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + + # u + for i in 1:J*N + out[i]=d0[i]+dot(FU,d1[i])+dot(FU,d2[i]*FU)-u[i] + end + # e + out[J*N+1]=-e+2*pi*rho* + ((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/ + (1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU))) + + return out +end +# DXi +@everywhere function simpleq_hardcore_DXi(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + out=Array{Float64,2}(undef,J*N+1,J*N+1) + FU=chebyshev(u,taus,weights,P,N,J,4) + + #D's + d0=D0(e,rho,r,weights,N,J) + d1=D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + d2=D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + dsed0=dseD0(e,rho,r,weights,N,J) + dsed1=dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + dsed2=dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + + # denominator of e + denom=1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU)) + + for zetapp in 0:J-1 + for n in 1:N + one=zeros(Int64,J*N) + one[zetapp*N+n]=1 + Fone=chebyshev(one,taus,weights,P,N,J,4) + + for i in 1:J*N + # du/du + out[i,zetapp*N+n]=dot(Fone,d1[i])+2*dot(FU,d2[i]*Fone)-(zetapp*N+n==i ? 1 : 0) + # du/de + out[i,J*N+1]=(dsed0[i]+dot(FU,dsed1[i])+dot(FU,dsed2[i]*FU))/(2*sqrt(abs(e)))*(e>=0 ? 1 : -1) + end + # de/du + out[J*N+1,zetapp*N+n]=2*pi*rho* + (-dot(Fone,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-2*dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*Fone))/denom- + 2*pi*rho* + ((1+2*sqrt(abs(e)))-gamma0(e,rho,weights)-dot(FU,gamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,gamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))* + rho^2*(dot(Fone,gammabar1(taus,T,weights,P,J,4))+2*dot(FU,gammabar2(taus,T,weights,P,J,4)*Fone))/denom^2 + end + end + #de/de + out[J*N+1,J*N+1]=-1+2*pi*rho* + (2-dsedgamma0(e,rho,weights)-dot(FU,dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))-dot(FU,dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/denom/ + (2*sqrt(abs(e)))*(e>=0 ? 1 : -1) + + return out +end + +# dXi/dmu +@everywhere function simpleq_hardcore_dXidmu(u,e,rho,r,taus,T,weights,weights_gL,P,N,J) + out=Array{Float64}(undef,J*N+1) + FU=chebyshev(u,taus,weights,P,N,J,4) + + #D's + dsmud0=dsmuD0(e,rho,r,weights,N,J) + dsmud1=dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + dsmud2=dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,4) + + # u + for i in 1:J*N + out[i]=(dsmud0[i]+dot(FU,dsmud1[i])+dot(FU,dsmud2[i]*FU))/(4*sqrt(abs(e))) + end + # e + out[J*N+1]=2*pi*rho*(2/3+1/(2*sqrt(abs(e)))- + (dsmudgamma0(e,rho,weights)+dot(FU,dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,4))+dot(FU,dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,4)*FU))/(4*sqrt(abs(e))) + )/(1-8/3*pi*rho+rho^2*(gammabar0(weights)+dot(FU,gammabar1(taus,T,weights,P,J,4))+dot(FU,gammabar2(taus,T,weights,P,J,4)*FU))) + + return out +end + +# B's +@everywhere function B0(r) + return pi/12*(r-1)^2*(r+5) +end +@everywhere function B1(r,zeta,n,taus,T,weights,nu) + return (taus[zeta+1]>=(2-r)/r || taus[zeta+2]<=-r/(r+2) ? 0 : + 8*pi/(r+1)*integrate_legendre(tau-> + (1-(r-(1-tau)/(1+tau))^2)/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1])),max(taus[zeta+1],-r/(r+2)) + ,min(taus[zeta+2],(2-r)/r),weights) + ) +end +@everywhere function B2(r,zeta,n,zetap,m,taus,T,weights,nu) + return 32*pi/(r+1)*integrate_legendre(tau-> + 1/(1+tau)^(3-nu)*T[n+1]((2*tau-(taus[zeta+1]+taus[zeta+2]))/(taus[zeta+2]-taus[zeta+1]))* + (taus[zetap+1]>=alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau) || taus[zetap+2]<=alpham(1+r,tau) ? 0 : + integrate_legendre(s-> + 1/(1+s)^(3-nu)*T[m+1]((2*s-(taus[zetap+1]+taus[zetap+2]))/(taus[zetap+2]-taus[zetap+1])),max(taus[zetap+1] + ,alpham(1+r,tau)),min(taus[zetap+2],alphap(abs(r-(1-tau)/(1+tau))-2*tau/(1+tau),tau)),weights) + ) + ,taus[zeta+1],taus[zeta+2],weights) +end + +# D's +@everywhere function D0(e,rho,r,weights,N,J) + out=Array{Float64}(undef,J*N) + for i in 1:J*N + out[i]=exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ + rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> + (s+1)*B0(s)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 + ,0,min(1,r[i]),weights)+ + (r[i]>=1 ? 0 : + rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_legendre(s-> + (s+r[i]+1)*B0(s+r[i])*exp(-2*sqrt(abs(e))*s) + ,0,1-r[i],weights) + ) + end + return out +end +@everywhere function D1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64}(undef,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> + (s+1)*B1(s,zeta,n,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 + ,0,r[i],weights)+ + rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s-> + (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + + return out +end +@everywhere function D2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64,2}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho*sqrt(abs(e))/(r[i]+1)*integrate_legendre(s-> + (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2 + ,0,r[i],weights)+ + rho*sqrt(abs(e))/(2*(r[i]+1))*(1-exp(-4*sqrt(abs(e))*r[i]))*integrate_laguerre(s-> + (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + end + end + return out +end + +# dD/d sqrt(abs(e))'s +@everywhere function dseD0(e,rho,r,weights,N,J) + out=Array{Float64}(undef,J*N) + for i in 1:J*N + out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ + rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B0(s)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,min(1,r[i]),weights)+ + (r[i]>=1 ? 0 : + rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights) + ) + end + return out +end +@everywhere function dseD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64}(undef,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + return out +end +@everywhere function dseD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64,2}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + end + end + return out +end + +# dD/d sqrt(abs(e+mu/2))'s +@everywhere function dsmuD0(e,rho,r,weights,N,J) + out=Array{Float64}(undef,J*N) + for i in 1:J*N + out[i]=-2*r[i]*exp(-2*sqrt(abs(e))*r[i])/(r[i]+1)+ + rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B0(s)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,min(1,r[i]),weights)+ + (r[i]>=1 ? 0 : + rho/(2*(r[i]+1))*integrate_legendre(s->(s+r[i]+1)*B0(s+r[i])*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))*exp(-2*sqrt(abs(e))*s)+4*sqrt(abs(e))*r[i]*exp(-2*sqrt(abs(e))*(2*r[i]+s))),0,1-r[i],weights) + ) + end + return out +end +@everywhere function dsmuD1(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64}(undef,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + out[i][zeta*(P+1)+n+1]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B1(s,zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B1(s+r[i],zeta,n,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + return out +end +@everywhere function dsmuD2(e,rho,r,taus,T,weights,weights_gL,P,N,J,nu) + out=Array{Array{Float64,2}}(undef,J*N) + for i in 1:J*N + out[i]=Array{Float64,2}(undef,(P+1)*J,(P+1)*J) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[i][zeta*(P+1)+n+1,zetap*(P+1)+m+1]=rho/(r[i]+1)*integrate_legendre(s-> + (s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*r[i])*(exp(-2*sqrt(abs(e))*(r[i]-s))-exp(-2*sqrt(abs(e))*(r[i]+s)))/2+2*sqrt(abs(e))*s*(exp(-2*sqrt(abs(e))*(r[i]-s))+exp(-2*sqrt(abs(e))*(r[i]+s)))/2) + ,0,r[i],weights)+ + rho/(2*(r[i]+1))*integrate_laguerre(s-> + (s+r[i]+1)*B2(s+r[i],zeta,n,zetap,m,taus,T,weights,nu)*((-1-2*sqrt(abs(e))*s)*(1-exp(-4*sqrt(abs(e))*r[i]))+4*sqrt(abs(e))*r[i]*exp(-4*sqrt(abs(e))*r[i])) + ,2*sqrt(abs(e)),weights_gL) + end + end + end + end + end + return out +end + +# gamma's +@everywhere function gamma0(e,rho,weights) + return 2*rho*e*integrate_legendre(s->(s+1)*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights) +end +@everywhere function gamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=2*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + return out +end +@everywhere function gamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=2*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + end + end + return out +end + +# dgamma/d sqrt(abs(e))'s +@everywhere function dsedgamma0(e,rho,weights) + return 4*rho*sqrt(abs(e))*integrate_legendre(s->(s+1)*B0(s)*(1-sqrt(abs(e))*s)*exp(-2*sqrt(abs(e))*s),0,1,weights) +end +@everywhere function dsedgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=4*rho*e*integrate_laguerre(s->(s+1)*B1(s,zeta,n,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL) + end + end + return out +end +@everywhere function dsedgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*rho*e*integrate_laguerre(s->(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu)*(1-sqrt(abs(e))*s),2*sqrt(abs(e)),weights_gL) + end + end + end + end + return out +end + +# dgamma/d sqrt(e+mu/2)'s +@everywhere function dsmudgamma0(e,rho,weights) + return -4*rho*e*integrate_legendre(s->(s+1)*s*B0(s)*exp(-2*sqrt(abs(e))*s),0,1,weights) +end +@everywhere function dsmudgamma1(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B1(s,zeta,n,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + return out +end +@everywhere function dsmudgamma2(e,rho,taus,T,weights,weights_gL,P,J,nu) + out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=-4*rho*e*integrate_laguerre(s->s*(s+1)*B2(s,zeta,n,zetap,m,taus,T,weights,nu),2*sqrt(abs(e)),weights_gL) + end + end + end + end + return out +end + +# \bar gamma's +@everywhere function gammabar0(weights) + return 4*pi*integrate_legendre(s->s^2*B0(s-1),0,1,weights) +end +@everywhere function gammabar1(taus,T,weights,P,J,nu) + out=Array{Float64}(undef,J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + out[zeta*(P+1)+n+1]=4*pi*integrate_legendre(s->s^2*B1(s-1,zeta,n,taus,T,weights,nu),0,1,weights) + end + end + return out +end +@everywhere function gammabar2(taus,T,weights,P,J,nu) + out=Array{Float64,2}(undef,J*(P+1),J*(P+1)) + for zeta in 0:J-1 + for n in 0:P + for zetap in 0:J-1 + for m in 0:P + out[zeta*(P+1)+n+1,zetap*(P+1)+m+1]=4*pi*integrate_legendre(s->s^2*B2(s-1,zeta,n,zetap,m,taus,T,weights,nu),0,1,weights) + end + end + end + end + return out +end diff --git a/src/simpleq-iteration.jl b/src/simpleq-iteration.jl new file mode 100644 index 0000000..98977b8 --- /dev/null +++ b/src/simpleq-iteration.jl @@ -0,0 +1,94 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# compute rho(e) using the iteration +function simpleq_iteration_rho_e(es,order,v,maxiter) + for j in 1:length(es) + (u,rho)=simpleq_iteration_hatun(es[j],order,v,maxiter) + @printf("% .15e % .15e\n",es[j],real(rho[maxiter+1])) + end +end + +# compute u(x) using the iteration and print at every step +function simpleq_iteration_ux(order,e,v,maxiter,xmin,xmax,nx) + (u,rho)=simpleq_iteration_hatun(e,order,v,maxiter) + + weights=gausslegendre(order) + for i in 1:nx + x=xmin+(xmax-xmin)*i/nx + @printf("% .15e ",x) + for n in 2:maxiter+1 + @printf("% .15e ",real(easyeq_u_x(x,u[n],weights))) + end + print('\n') + end +end + + +# \hat u_n +function simpleq_iteration_hatun(e,order,v,maxiter) + # gauss legendre weights + weights=gausslegendre(order) + + # initialize V and Eta + (V,V0)=easyeq_init_v(weights,v) + (Eta,Eta0)=easyeq_init_H(weights,v) + + # init u and rho + u=Array{Array{Float64}}(undef,maxiter+1) + u[1]=zeros(Float64,order) + rho=zeros(Float64,maxiter+1) + + # iterate + for n in 1:maxiter + u[n+1]=simpleq_iteration_A(e,weights,Eta)\simpleq_iteration_b(u[n],e,rho[n],V) + rho[n+1]=simpleq_iteration_rhon(u[n+1],e,weights,V0,Eta0) + end + + return (u,rho) +end + +# A +function simpleq_iteration_A(e,weights,Eta) + N=length(weights[1]) + out=zeros(Float64,N,N) + for i in 1:N + k=(1-weights[1][i])/(1+weights[1][i]) + out[i,i]=k^2+4*e + for j in 1:N + y=(weights[1][j]+1)/2 + out[i,j]+=weights[2][j]*(1-y)*Eta[i][j]/(2*(2*pi)^3*y^3) + end + end + return out +end + +# b +function simpleq_iteration_b(u,e,rho,V) + out=zeros(Float64,length(V)) + for i in 1:length(V) + out[i]=V[i]+2*e*rho*u[i]^2 + end + return out +end + +# rho_n +function simpleq_iteration_rhon(u,e,weights,V0,Eta0) + S=V0 + for i in 1:length(weights[1]) + y=(weights[1][i]+1)/2 + S+=-weights[2][i]*(1-y)*u[i]*Eta0[i]/(2*(2*pi)^3*y^3) + end + return 2*e/S +end diff --git a/src/tools.jl b/src/tools.jl new file mode 100644 index 0000000..0d3dc7f --- /dev/null +++ b/src/tools.jl @@ -0,0 +1,49 @@ +## Copyright 2021 Ian Jauslin +## +## Licensed under the Apache License, Version 2.0 (the "License"); +## you may not use this file except in compliance with the License. +## You may obtain a copy of the License at +## +## http://www.apache.org/licenses/LICENSE-2.0 +## +## Unless required by applicable law or agreed to in writing, software +## distributed under the License is distributed on an "AS IS" BASIS, +## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +## See the License for the specific language governing permissions and +## limitations under the License. + +# \Phi(x)=2*(1-sqrt(1-x))/x +@everywhere function Phi(x) + if abs(x)>1e-5 + return 2*(1-sqrt(abs(1-x)))/x + else + return 1+x/4+x^2/8+5*x^3/64+7*x^4/128+21*x^5/512 + end +end +# \partial\Phi +@everywhere function dPhi(x) + #if abs(x-1)<1e-5 + # @printf(stderr,"warning: dPhi is singular at 1, and evaluating it at (% .8e+i% .8e)\n",real(x),imag(x)) + #end + if abs(x)>1e-5 + return 1/(sqrt(abs(1-x))*x)*(1-x>=0 ? 1 : -1)-2*(1-sqrt(abs(1-x)))/x^2 + else + return 1/4+x/4+15*x^2/64+7*x^3/32+105*x^4/512+99*x^5/512 + end +end + +# apply Phi to every element of a vector +@everywhere function dotPhi(v) + out=zeros(Float64,length(v)) + for i in 1:length(v) + out[i]=Phi(v[i]) + end + return out +end +@everywhere function dotdPhi(v) + out=zeros(Float64,length(v)) + for i in 1:length(v) + out[i]=dPhi(v[i]) + end + return out +end |