diff options
Diffstat (limited to 'src/main.jl')
| -rw-r--r-- | src/main.jl | 168 |
1 files changed, 133 insertions, 35 deletions
diff --git a/src/main.jl b/src/main.jl index 382fe6b..28fc2be 100644 --- a/src/main.jl +++ b/src/main.jl @@ -1,4 +1,4 @@ -## Copyright 2021 Ian Jauslin +## Copyright 2021-2023 Ian Jauslin ## ## Licensed under the Apache License, Version 2.0 (the "License"); ## you may not use this file except in compliance with the License. @@ -22,6 +22,8 @@ include("chebyshev.jl") include("integration.jl") include("interpolation.jl") include("tools.jl") +include("multithread.jl") +include("optimization.jl") include("potentials.jl") include("print.jl") include("easyeq.jl") @@ -36,9 +38,6 @@ function main() rho=1e-6 e=1e-4 - # incrementally initialize Newton algorithm - nlrho_init=1 - # potential v=k->v_exp(k,1.) a0=a0_exp(1.) @@ -49,19 +48,29 @@ function main() v_param_e=1. # plot range when plotting in rho - minlrho=-6 - maxlrho=2 + # linear + minrho=1e-6 + maxrho=1e2 + nrho=0 + # logarithmic + minlrho=-6. + maxlrho=2. nlrho=100 - rhos=Array{Float64}(undef,0) + # list + rhos=Array{Float64,1}(undef,0) # plot range when plotting in e - minle=-6 - maxle=2 + minle=-6. + maxle=2. nle=100 - es=Array{Float64}(undef,0) + es=Array{Float64,1}(undef,0) # plot range when plotting in x - xmin=0 - xmax=100 + xmin=0. + xmax=100. nx=100 + # plot range when plotting in k + kmin=0. + kmax=10. + nk=100 # cutoffs tolerance=1e-11 @@ -77,8 +86,8 @@ function main() 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 + minlrho_init=-6. + nlrho_init=0 # Hann window for Fourier transforms windowL=1e3 @@ -98,6 +107,18 @@ function main() anyeq_compleq_approx=Anyeq_approx(1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.) anyeq_approx=anyeq_bigeq_approx + # numerical approximations of derivatives + dx=1e-7 + dk=1e-7 + + # initial guess for 2pt_max + x0=1. + k0=1. + # maximum step in 2pt_max + maxstep=Inf + # tolerance for max search + tolerance_max=Inf + # read cli arguments (params,potential,method,savefile,command)=read_args(ARGS) @@ -109,8 +130,8 @@ function main() print(stderr,"error: could not read parameter '",param,"'.\n") exit(-1) end - lhs=terms[1] - rhs=terms[2] + lhs=string(terms[1]) + rhs=string(terms[2]) if lhs=="rho" rho=parse(Float64,rhs) elseif lhs=="minlrho_init" @@ -133,6 +154,12 @@ function main() maxlrho=parse(Float64,rhs) elseif lhs=="nlrho" nlrho=parse(Int64,rhs) + elseif lhs=="minrho" + minrho=parse(Float64,rhs) + elseif lhs=="maxrho" + maxrho=parse(Float64,rhs) + elseif lhs=="nrho" + nrho=parse(Int64,rhs) elseif lhs=="es" es=parse_list(rhs) elseif lhs=="minle" @@ -147,6 +174,12 @@ function main() xmax=parse(Float64,rhs) elseif lhs=="nx" nx=parse(Int64,rhs) + elseif lhs=="kmin" + kmin=parse(Float64,rhs) + elseif lhs=="kmax" + kmax=parse(Float64,rhs) + elseif lhs=="nk" + nk=parse(Int64,rhs) elseif lhs=="P" P=parse(Int64,rhs) elseif lhs=="N" @@ -155,6 +188,18 @@ function main() J=parse(Int64,rhs) elseif lhs=="window_L" windowL=parse(Float64,rhs) + elseif lhs=="dx" + dx=parse(Float64,rhs) + elseif lhs=="x0" + x0=parse(Float64,rhs) + elseif lhs=="dk" + dk=parse(Float64,rhs) + elseif lhs=="k0" + k0=parse(Float64,rhs) + elseif lhs=="maxstep" + maxstep=parse(Float64,rhs) + elseif lhs=="tolerance_max" + tolerance_max=parse(Float64,rhs) elseif lhs=="aK" anyeq_approx.aK=parse(Float64,rhs) elseif lhs=="bK" @@ -242,37 +287,48 @@ function main() ## 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)) + # linear only if nrho is specified + if nrho>0 + rhos=Array{Float64,1}(undef,nrho) + for j in 0:nrho-1 + rhos[j+1]=(nrho==1 ? minrho : minrho+(maxrho-minrho)/(nrho-1)*j) + end + else + rhos=Array{Float64,1}(undef,nlrho) + for j in 0:nlrho-1 + rhos[j+1]=(nlrho==1 ? 10^minlrho : 10^(minlrho+(maxlrho-minlrho)/(nlrho-1)*j)) + end end end + # es if length(es)==0 - es=Array{Float64}(undef,nle) + es=Array{Float64,1}(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) + taus=Array{Float64,1}(undef,J+1) for j in 0:J taus[j+1]=-1+2*j/J end + # tolerance_max + if tolerance_max==Inf + tolerance_max=tolerance + end + ## run command - if method=="easyeq" + if command=="scattering_length" + @printf("% .15e\n",a0) + elseif 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) + easyeq_energy_rho(rhos,minlrho_init,nlrho_init,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) @@ -285,7 +341,28 @@ function main() 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) + easyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,order,a0,v,maxiter,tolerance,easyeq_approx) + # 2pt correlation + elseif command=="2pt" + easyeq_2pt(xmin,xmax,nx,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx) + # max of 2pt correlation + elseif command=="2pt_max" + easyeq_2pt_max(dx,x0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # max of 2pt correlation as a function of rho + elseif command=="2pt_max_rho" + easyeq_2pt_max_rho(rhos,dx,x0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # fourier transform of 2pt correlation + elseif command=="2pt_fourier" + easyeq_2pt_fourier(kmin,kmax,nk,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx) + # max of fourier transform of 2pt correlation + elseif command=="2pt_fourier_max" + easyeq_2pt_fourier_max(dk,k0,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # max of 2pt correlation as a function of rho + elseif command=="2pt_fourier_max_rho" + easyeq_2pt_fourier_max_rho(rhos,dk,k0,minlrho_init,nlrho_init,order,windowL,a0,v,maxstep,maxiter,tolerance,tolerance_max,easyeq_approx) + # momentum distribution + elseif command=="momentum_distribution" + easyeq_momentum_distribution(kmin,kmax,minlrho_init,nlrho_init,order,windowL,rho,a0,v,maxiter,tolerance,easyeq_approx) else print(stderr,"unrecognized command '",command,"'.\n") exit(-1) @@ -316,7 +393,7 @@ function main() 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) + anyeq_energy_rho(rhos,minlrho_init,nlrho_init,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" @@ -331,12 +408,29 @@ function main() 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) + anyeq_condensate_fraction_rho(rhos,minlrho_init,nlrho_init,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) + anyeq_momentum_distribution(kmin,kmax,rho,minlrho_init,nlrho_init,taus,P,N,J,windowL,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) + elseif command=="2pt_max" + anyeq_2pt_correlation_max(rho,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="2pt_max_rho" + anyeq_2pt_correlation_max_rho(rhos,minlrho_init,nlrho_init,dx,x0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="2pt_fourier" + anyeq_2pt_correlation_fourier(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,kmin,kmax,nk,anyeq_approx,savefile) + elseif command=="2pt_fourier_test" + anyeq_2pt_correlation_fourier_test(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmax,kmin,kmax,nk,anyeq_approx,savefile) + elseif command=="2pt_fourier_max" + anyeq_2pt_correlation_fourier_max(rho,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="2pt_fourier_max_rho" + anyeq_2pt_correlation_fourier_max_rho(rhos,minlrho_init,nlrho_init,dk,k0,maxstep,taus,P,N,J,windowL,a0,v,maxiter,tolerance,tolerance_max,anyeq_approx,savefile) + elseif command=="uncondensed_2pt" + anyeq_uncondensed_2pt_correlation(minlrho_init,nlrho_init,taus,P,N,J,windowL,rho,a0,v,maxiter,tolerance,xmin,xmax,nx,anyeq_approx,savefile) + # compressibility + elseif command=="compressibility_rho" + anyeq_compressibility_rho(rhos,minlrho_init,nlrho_init,taus,P,N,J,a0,v,maxiter,tolerance,anyeq_approx,savefile) else print(stderr,"unrecognized command: '",command,"'\n") exit(-1) @@ -360,9 +454,11 @@ function main() end # parse a comma separated list as an array of Float64 -function parse_list(str) +function parse_list( + str::String +) elems=split(str,",") - out=Array{Float64}(undef,length(elems)) + out=Array{Float64,1}(undef,length(elems)) for i in 1:length(elems) out[i]=parse(Float64,elems[i]) end @@ -370,7 +466,9 @@ function parse_list(str) end # read cli arguments -function read_args(ARGS) +function read_args( + ARGS +) # flag flag="" |