From d096cbb1007449c2847b93279ea7a476eed135c3 Mon Sep 17 00:00:00 2001 From: Ian Jauslin Date: Thu, 6 Apr 2023 11:29:01 -0400 Subject: remove julia --- julia/driving.jl | 17 ---- julia/main.jl | 150 ---------------------------------- julia/navier_stokes.jl | 215 ------------------------------------------------- julia/print.jl | 15 ---- 4 files changed, 397 deletions(-) delete mode 100644 julia/driving.jl delete mode 100644 julia/main.jl delete mode 100644 julia/navier_stokes.jl delete mode 100644 julia/print.jl diff --git a/julia/driving.jl b/julia/driving.jl deleted file mode 100644 index 657f252..0000000 --- a/julia/driving.jl +++ /dev/null @@ -1,17 +0,0 @@ -# driving forces (g stands for k^\perp.g/(2i\pi|k|) - -# a simple test driving force -function g_test( - kx::Int64, - ky::Int64 -) - if(kx==2 && ky==-1) - return 0.5+sqrt(3)/2*1im - elseif (kx==-2 && ky==1) - return 0.5-sqrt(3)/2*1im - else - return 0. - end -end - - diff --git a/julia/main.jl b/julia/main.jl deleted file mode 100644 index 33444f3..0000000 --- a/julia/main.jl +++ /dev/null @@ -1,150 +0,0 @@ -using Printf -using FFTW -include("navier_stokes.jl") -include("driving.jl") -include("print.jl") - -function main() - ## defaults - delta=0.0001220703125 - nsteps=10000000 - nu=0.00048828125 - print_freq=1000 - - K1=16 - K2=K1 - N1=4*K1+1 - N2=4*K2+1 - - # read cli arguments - (params,driving,command)=read_args(ARGS) - - # whether N was set in parameters - setN1=false - setN2=false - - # 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=string(terms[1]) - rhs=string(terms[2]) - if lhs=="K1" - K1=parse(Int64,rhs) - elseif lhs=="K2" - K2=parse(Int64,rhs) - elseif lhs=="K" - K1=parse(Int64,rhs) - K2=K1 - elseif lhs=="N1" - N1=parse(Int64,rhs) - setN1=true - elseif lhs=="N2" - N2=parse(Int64,rhs) - setN2=true - elseif lhs=="N" - N1=parse(Int64,rhs) - N2=N1 - setN1=true - setN2=true - elseif lhs=="nsteps" - nsteps=parse(Int64,rhs) - elseif lhs=="nu" - nu=parse(Float64,rhs) - elseif lhs=="delta" - delta=parse(Float64,rhs) - elseif lhs=="print_freq" - print_freq=parse(Int64,rhs) - else - print(stderr,"unrecognized parameter '",lhs,"'.\n") - exit(-1) - end - end - end - - ## read driving - if driving=="test" - g=(kx,ky)->g_test(kx,ky) - end - - ## set parameters - if ! setN1 - N1=4*K1+1 - end - if ! setN2 - N2=4*K2+1 - end - - ## run command - if command=="uk" - navier_stokes_uk(nsteps,nu,g,N1,N2,K1,K2,delta,print_freq) - elseif command=="alpha" - navier_stokes_alpha(nsteps,nu,g,N1,N2,K1,K2,delta,print_freq) - else - print(stderr,"unrecognized command '",command,"'.\n") - exit(-1) - end -end - -# read cli arguments -function read_args(ARGS) - # flag - flag="" - - # output strings - command="" - params="" - # default driving force - driving="test" - - # 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=='g' - # raise flag - flag="driving" - else - print_usage() - exit(-1) - end - end - # arguments that do not start with a dash - else - if flag=="params" - params=arg - elseif flag=="driving" - driving=arg - else - command=arg - end - # reset flag - flag="" - end - end - - if command=="" - print_usage() - exit(-1) - end - - return (params,driving,command) -end - -# usage -function print_usage() - print(stderr,"usage: nstrophy [-p params] [-g driving] \n") -end - -main() diff --git a/julia/navier_stokes.jl b/julia/navier_stokes.jl deleted file mode 100644 index a76e918..0000000 --- a/julia/navier_stokes.jl +++ /dev/null @@ -1,215 +0,0 @@ -# solve Navier-Stokes -function navier_stokes_uk( - nsteps::Int64, - nu::Float64, - g::Function, - N1::Int64, - N2::Int64, - K1::Int64, - K2::Int64, - delta::Float64, - print_freq::Int64 -) - # init - u=ns_init(K1,K2) - - for i in 1:nsteps - u=ns_RK4(u,nu,g,N1,N2,K1,K2,delta) - - if i % print_freq==0 - @printf("% .15e ",delta*i) - print_u_inline(u,N1,N2,K1,K2) - - # print to stderr - @printf(stderr,"% .8e ",delta*i) - for k1 in -K1:K1 - for k2 in -K2:K2 - index=indices_from_momentum(k1,2*K1+1)*(2*K2+1)+indices_from_momentum(k2,2*K2+1)+1 - if (k1^2+k2^2<=1) - @printf(stderr,"% .8e % .8e ",real(u[index]),imag(u[index])) - end - end - end - @printf(stderr,"\n") - end - end - -end - -# compute the enstrophy -function navier_stokes_alpha( - nsteps::Int64, - nu::Float64, - g::Function, - N1::Int64, - N2::Int64, - K1::Int64, - K2::Int64, - delta::Float64, - print_freq::Int64 -) - # init - u=ns_init(K1,K2) - - for i in 1:nsteps - u=ns_RK4(u,nu,g,N1,N2,K1,K2,delta) - if i % print_freq==0 - alpha=ns_alpha(u,g,K1,K2) - @printf("% .15e % .15e % .15e\n",i*delta,real(alpha),imag(alpha)) - @printf(stderr,"% .15e % .15e % .15e\n",i*delta,real(alpha),imag(alpha)) - end - end -end - -# initialize u -function ns_init( - K1::Int64, - K2::Int64 -) - u=zeros(Complex{Float64},(2*K1+1)*(2*K2+1)) - for k1 in -K1:K1 - for k2 in -K2:K2 - index=indices_from_momentum(k1,2*K1+1)*(2*K2+1)+indices_from_momentum(k2,2*K2+1)+1 - u[index]=exp(-sqrt(k1^2+k2^2)) - end - end - - return u -end - -# Navier-Stokes equation (right side) -function ns_rhs( - u::Array{Complex{Float64},1}, - nu::Float64, - g::Function, - N1::Int64, - N2::Int64, - K1::Int64, - K2::Int64 -) - # compute the T term - out=ns_T(u,N1,N2,K1,K2) - - # add the rest of the terms - for k1 in -K1:K1 - for k2 in -K2:K2 - index=indices_from_momentum(k1,2*K1+1)*(2*K2+1)+indices_from_momentum(k2,2*K2+1)+1 - out[index]=-(4*pi^2*nu*(k1^2+k2^2))*u[index]+g(k1,k2)+(k1^2+k2^2>0 ? 4*pi^2/sqrt(k1^2+k2^2)*out[index] : 0) - end - end - - return out -end - -# compute k\in K from n\in N -function momentum_from_indices( - i::Int64, - K::Int64, - N::Int64, -) - if i<=K - k=i - elseif i>=N-K - k=i-N - else - k=i - end - - return(k) -end -# inverse -function indices_from_momentum( - k::Int64, - N::Int64 -) - if k>=0 - return k - else - return N+k - end -end - -# T term in Navier-Stokes -function ns_T( - u::Array{Complex{Float64},1}, - N1::Int64, - N2::Int64, - K1::Int64, - K2::Int64 -) - Fx1=zeros(Complex{Float64},N1,N2) - Fy1=zeros(Complex{Float64},N1,N2) - Fx2=zeros(Complex{Float64},N1,N2) - Fy2=zeros(Complex{Float64},N1,N2) - - for k1 in -K1:K1 - for k2 in -K2:K2 - i1=indices_from_momentum(k1,N1)+1 - i2=indices_from_momentum(k2,N2)+1 - index=indices_from_momentum(k1,2*K1+1)*(2*K2+1)+indices_from_momentum(k2,2*K2+1)+1 - # no need to divide by N: ifft does the division already - Fx1[i1,i2]=(k1==0 ? 0. : k1/sqrt(k1^2+k2^2)*u[index]) - Fx2[i1,i2]=(k1==0 ? 0. : k1*sqrt(k1^2+k2^2)*u[index]) - Fy1[i1,i2]=(k2==0 ? 0. : k2/sqrt(k1^2+k2^2)*u[index]) - Fy2[i1,i2]=(k2==0 ? 0. : k2*sqrt(k1^2+k2^2)*u[index]) - end - end - - fft!(Fx1) - fft!(Fx2) - fft!(Fy1) - fft!(Fy2) - Fx1=Fx1.*Fy2.-Fx2.*Fy1 - ifft!(Fx1) - - out=zeros(Complex{Float64},(2*K1+1)*(2*K2+1)) - for k1 in -K1:K1 - for k2 in -K2:K2 - i1=indices_from_momentum(k1,N1)+1 - i2=indices_from_momentum(k2,N2)+1 - index=indices_from_momentum(k1,2*K1+1)*(2*K2+1)+indices_from_momentum(k2,2*K2+1)+1 - out[index]=Fx1[i1,i2] - end - end - - return out -end - -# Runge-Kutta step to solve Navier-Stokes -function ns_RK4( - u::Array{Complex{Float64},1}, - nu::Float64, - g::Function, - N1::Int64, - N2::Int64, - K1::Int64, - K2::Int64, - delta::Float64 -) - k1=delta*ns_rhs(u,nu,g,N1,N2,K1,K2) - k2=delta*ns_rhs(u.+(k1./2),nu,g,N1,N2,K1,K2) - k3=delta*ns_rhs(u.+(k2./2),nu,g,N1,N2,K1,K2) - k4=delta*ns_rhs(u.+k3,nu,g,N1,N2,K1,K2) - - return u.+(k1+2 .*k2+2 .*k3.+k4)./6 -end - -# enstrophy -function ns_alpha( - u::Array{Complex{Float64},1}, - g::Function, - K1::Int64, - K2::Int64 -) - denom=0. *1im - num=0. *1im - for k1 in -K1:K1 - for k2 in -K2:K2 - index=indices_from_momentum(k1,2*K1+1)*(2*K2+1)+indices_from_momentum(k2,2*K2+1)+1 - denom+=(k1^2+k2^2)^2*u[index]*conj(u[index])*(1. +(k2!=0 ? k1^2/k2^2 : 0.)) - num+=(k1^2+k2^2)*u[index]*conj(g(k1,k2))*(1. +(k2!=0 ? k1^2/k2^2 : 0.)) - end - end - - return num/denom -end diff --git a/julia/print.jl b/julia/print.jl deleted file mode 100644 index 61be40d..0000000 --- a/julia/print.jl +++ /dev/null @@ -1,15 +0,0 @@ -function print_u_inline( - u::Array{Complex{Float64},1}, - N1::Int64, - N2::Int64, - K1::Int64, - K2::Int64 -) - for k1 in -K1:K1 - for k2 in -K2:K2 - index=indices_from_momentum(k1,2*K1+1)*(K2+1)+indices_from_momentum(k2,2*K2+1)+1 - @printf("% .15e % .15e ",real(u[index]),imag(u[index])) - end - end - print("\n") -end -- cgit v1.2.3-54-g00ecf