diff options
Diffstat (limited to 'src/anyeq.jl')
| -rw-r--r-- | src/anyeq.jl | 949 |
1 files changed, 949 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 |