Skip to content

Commit

Permalink
changes variable names in TwoDTurb from U, V --> u, v
Browse files Browse the repository at this point in the history
  • Loading branch information
navidcy committed Dec 29, 2018
1 parent 8b4dc93 commit 8280190
Show file tree
Hide file tree
Showing 5 changed files with 44 additions and 50 deletions.
2 changes: 1 addition & 1 deletion examples/twodturb_randomdecay.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ while cl.step < ntot
end

# Plot the radial energy spectrum
E = @. 0.5*(vs.U^2 + vs.V^2) # energy density
E = @. 0.5*(vs.u^2 + vs.v^2) # energy density
Eh = rfft(E)
kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1)

Expand Down
12 changes: 6 additions & 6 deletions examples/twodturb_stochasticforcing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ function makeplot(prob, diags)

t = round(mu*cl.t, digits=2)
sca(axs[1]); cla()
pcolormesh(x, y, prob.vars.q)
pcolormesh(x, y, v.q)
xlabel(L"$x$")
ylabel(L"$y$")
title("\$\\nabla^2\\psi(x,y,\\mu t= $t )\$")
Expand All @@ -74,15 +74,15 @@ function makeplot(prob, diags)

# If the Ito interpretation was used for the work
# then we need to add the drift term
# total = W[ii2]+σ - D[ii] - R[ii] # Ito
# total = W[ii2]+ε - D[ii] - R[ii] # Ito
total = W[ii2] - D[ii] - R[ii] # Stratonovich
residual = dEdt - total

# If the Ito interpretation was used for the work
# then we need to add the drift term: I[ii2] + σ
# then we need to add the drift term: I[ii2] + ε
plot(mu*E.t[ii], W[ii2], label=L"work ($W$)") # Ito
# plot(mu*E.t[ii], W[ii2] , label=L"work ($W$)") # Stratonovich
plot(mu*E.t[ii], σ .+ 0*E.t[ii], "--", label=L"ensemble mean work ($\langle W\rangle $)")
plot(mu*E.t[ii], ε .+ 0*E.t[ii], "--", label=L"ensemble mean work ($\langle W\rangle $)")
# plot(mu*E.t[ii], -D[ii], label="dissipation (\$D\$)")
plot(mu*E.t[ii], -R[ii], label=L"drag ($D=2\mu E$)")
plot(mu*E.t[ii], 0*E.t[ii], "k:", linewidth=0.5)
Expand Down Expand Up @@ -116,14 +116,14 @@ for i = 1:ns
TwoDTurb.updatevars!(prob)
# saveoutput(out)

cfl = cl.dt*maximum([maximum(v.U)/g.dx, maximum(v.V)/g.dy])
cfl = cl.dt*maximum([maximum(v.u)/g.dx, maximum(v.v)/g.dy])
res = makeplot(prob, diags)
pause(0.01)

@printf("step: %04d, t: %.1f, cfl: %.3f, time: %.2f s\n", cl.step, cl.t,
cfl, (time()-startwalltime)/60)

# savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), prob.step)
# savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
# savefig(savename, dpi=240)
end

Expand Down
60 changes: 30 additions & 30 deletions src/twodturb.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ using FourierFlows: getfieldspecs, structvarsexpr, parsevalsum, parsevalsum2

abstract type TwoDTurbVars <: AbstractVars end

const physicalvars = [:q, :U, :V]
const physicalvars = [:q, :u, :v]
const transformvars = [ Symbol(var, :h) for var in physicalvars ]
const forcedvars = [:Fh]
const stochforcedvars = [:prevsol]
Expand Down Expand Up @@ -120,9 +120,9 @@ eval(structvarsexpr(:StochasticForcedVars, stochforcedvarspecs; parent=:TwoDTurb
Returns the vars for unforced two-dimensional turbulence with grid g.
"""
function Vars(g; T=typeof(g.Lx))
@createarrays T (g.nx, g.ny) q U V
@createarrays Complex{T} (g.nkr, g.nl) sol qh Uh Vh
Vars(q, U, V, qh, Uh, Vh)
@createarrays T (g.nx, g.ny) q u v
@createarrays Complex{T} (g.nkr, g.nl) sol qh uh vh
Vars(q, u, v, qh, uh, vh)
end

"""
Expand Down Expand Up @@ -159,21 +159,21 @@ end
Calculates the advection term.
"""
function calcN_advection!(N, sol, t, cl, v, p, g)
@. v.Uh = im * g.l * g.invKrsq * sol
@. v.Vh = -im * g.kr * g.invKrsq * sol
@. v.uh = im * g.l * g.invKrsq * sol
@. v.vh = -im * g.kr * g.invKrsq * sol
@. v.qh = sol

ldiv!(v.U, g.rfftplan, v.Uh)
ldiv!(v.V, g.rfftplan, v.Vh)
ldiv!(v.u, g.rfftplan, v.uh)
ldiv!(v.v, g.rfftplan, v.vh)
ldiv!(v.q, g.rfftplan, v.qh)

@. v.U *= v.q # U*q
@. v.V *= v.q # V*q
@. v.u *= v.q # U*q
@. v.v *= v.q # V*q

mul!(v.Uh, g.rfftplan, v.U) # \hat{U*q}
mul!(v.Vh, g.rfftplan, v.V) # \hat{U*q}
mul!(v.uh, g.rfftplan, v.u) # \hat{U*q}
mul!(v.vh, g.rfftplan, v.v) # \hat{U*q}

@. N = -im*g.kr*v.Uh - im*g.l*v.Vh
@. N = -im*g.kr*v.uh - im*g.l*v.vh
nothing
end

Expand Down Expand Up @@ -213,11 +213,11 @@ Update the vars in v on the grid g with the solution in sol.
function updatevars!(prob)
v, g, sol = prob.vars, prob.grid, prob.sol
v.qh .= sol
@. v.Uh = im * g.l * g.invKrsq * sol
@. v.Vh = -im * g.kr * g.invKrsq * sol
@. v.uh = im * g.l * g.invKrsq * sol
@. v.vh = -im * g.kr * g.invKrsq * sol
ldiv!(v.q, g.rfftplan, deepcopy(v.qh))
ldiv!(v.U, g.rfftplan, deepcopy(v.Uh))
ldiv!(v.V, g.rfftplan, deepcopy(v.Vh))
ldiv!(v.u, g.rfftplan, deepcopy(v.uh))
ldiv!(v.v, g.rfftplan, deepcopy(v.vh))
nothing
end

Expand All @@ -243,8 +243,8 @@ solution `sol`.
"""
@inline function energy(prob)
sol, v, g = prob.sol, prob.vars, prob.grid
@. v.Uh = g.invKrsq * abs2(sol)
1/(2*g.Lx*g.Ly)*parsevalsum(v.Uh, g)
@. v.uh = g.invKrsq * abs2(sol)
1/(2*g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

"""
Expand All @@ -265,9 +265,9 @@ Returns the domain-averaged dissipation rate. nnu must be >= 1.
"""
@inline function dissipation(prob)
sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
@. v.Uh = g.Krsq^(p.nnu-1) * abs2(sol)
v.Uh[1, 1] = 0
p.nu/(g.Lx*g.Ly)*parsevalsum(v.Uh, g)
@. v.uh = g.Krsq^(p.nnu-1) * abs2(sol)
v.uh[1, 1] = 0
p.nu/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

"""
Expand All @@ -277,14 +277,14 @@ end
Returns the domain-averaged rate of work of energy by the forcing Fh.
"""
@inline function work(sol, v::ForcedVars, g)
@. v.Uh = g.invKrsq * sol * conj(v.Fh)
1/(g.Lx*g.Ly)*parsevalsum(v.Uh, g)
@. v.uh = g.invKrsq * sol * conj(v.Fh)
1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline function work(sol, v::StochasticForcedVars, g)
@. v.Uh = g.invKrsq * (v.prevsol + sol)/2.0 * conj(v.Fh) # Stratonovich
# @. v.Uh = g.invKrsq * v.prevsol * conj(v.Fh) # Ito
1/(g.Lx*g.Ly)*parsevalsum(v.Uh, g)
@. v.uh = g.invKrsq * (v.prevsol + sol)/2.0 * conj(v.Fh) # Stratonovich
# @. v.uh = g.invKrsq * v.prevsol * conj(v.Fh) # Ito
1/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

@inline work(prob) = work(prob.sol, prob.vars, prob.grid)
Expand All @@ -296,9 +296,9 @@ Returns the extraction of domain-averaged energy by drag/hypodrag mu.
"""
@inline function drag(prob)
sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
@. v.Uh = g.Krsq^(p.nmu-1) * abs2(sol)
v.Uh[1, 1] = 0
p.mu/(g.Lx*g.Ly)*parsevalsum(v.Uh, g)
@. v.uh = g.Krsq^(p.nmu-1) * abs2(sol)
v.uh[1, 1] = 0
p.mu/(g.Lx*g.Ly)*parsevalsum(v.uh, g)
end

end # module
4 changes: 2 additions & 2 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@ const rtol_niwqg = 1e-13 # tolerance for niwqg forcing tests
const rtol_twodturb = 1e-13 # tolerance for niwqg forcing tests

"Get the CFL number, assuming a uniform grid with `dx=dy`."
cfl(U, V, dt, dx) = maximum([maximum(abs.(U)), maximum(abs.(V))]*dt/dx)
cfl(prob) = cfl(prob.vars.U, prob.vars.V, prob.clock.dt, prob.grid.dx)
cfl(u, v, dt, dx) = maximum([maximum(abs.(u)), maximum(abs.(v))]*dt/dx)
cfl(prob) = cfl(prob.vars.u, prob.vars.v, prob.clock.dt, prob.grid.dx)

"Returns the energy in vertically Fourier mode 1 in the Boussinesq equations."
e1_fourier(u, v, p, m, N) = @. abs2(u) + abs2(v) + m^2*abs2(p)/N^2
Expand Down
16 changes: 5 additions & 11 deletions test/test_twodturb.jl
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ function test_twodturb_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7

# Forcing
kf, dkf = 12.0, 2.0
σ = 0.1
ε = 0.1

gr = TwoDGrid(n, L)
x, y = gridpoints(gr)
Expand All @@ -42,8 +42,8 @@ function test_twodturb_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7
@. force2k[gr.Krsq .< 2.0^2 ] = 0
@. force2k[gr.Krsq .> 20.0^2 ] = 0
@. force2k[Kr .< 2π/L] = 0
σ0 = parsevalsum(force2k.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly)
force2k .= σ/σ0 * force2k
ε0 = parsevalsum(force2k.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly)
force2k .= ε/ε0 * force2k

Random.seed!(1234)

Expand All @@ -70,7 +70,6 @@ function test_twodturb_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7
stepforward!(prob, diags, round(Int, nt))
TwoDTurb.updatevars!(prob)

cfl = cl.dt*maximum([maximum(v.V)/g.dx, maximum(v.U)/g.dy])
E, D, W, R = diags
t = round(mu*cl.t, digits=2)

Expand All @@ -82,15 +81,15 @@ function test_twodturb_stochasticforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7
# dEdt = W - D - R?
# If the Ito interpretation was used for the work
# then we need to add the drift term
# total = W[ii2]+σ - D[ii] - R[ii] # Ito
# total = W[ii2]+ε - D[ii] - R[ii] # Ito
total = W[ii2] - D[ii] - R[ii] # Stratonovich

residual = dEdt - total
isapprox(mean(abs.(residual)), 0, atol=1e-4)
end


function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, nmu=0, message=false)
function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1e-7, nnu=2, mu=1e-1, nmu=0)
n, L = 256, 2π
nu, nnu = 1e-7, 2
mu, nmu = 1e-1, 0
Expand Down Expand Up @@ -126,7 +125,6 @@ function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1
stepforward!(prob, diags, round(Int, nt))
TwoDTurb.updatevars!(prob)

cfl = cl.dt*maximum([maximum(v.V)/g.dx, maximum(v.U)/g.dy])
E, D, W, R = diags
t = round(mu*cl.t, digits=2)

Expand All @@ -140,10 +138,6 @@ function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1

residual = dEdt - total

if message
println("step: %04d, t: %.1f, cfl: %.3f, time: %.2f s\n", cl.step, cl.t, cfl, tc)
end

isapprox(mean(abs.(residual)), 0, atol=1e-8)
end

Expand Down

0 comments on commit 8280190

Please sign in to comment.