diff --git a/docs/src/modules/twodturb.md b/docs/src/modules/twodturb.md index f94a995d..b83105dd 100644 --- a/docs/src/modules/twodturb.md +++ b/docs/src/modules/twodturb.md @@ -9,11 +9,11 @@ This module solves two-dimensional incompressible turbulence. The flow is given through a streamfunction $\psi$ as $(u,\upsilon) = (-\partial_y\psi, \partial_x\psi)$. The dynamical variable used here is the component of the vorticity of the flow -normal to the plane of motion, $q=\partial_x \upsilon- \partial_y u = \nabla^2\psi$. +normal to the plane of motion, $\zeta=\partial_x \upsilon- \partial_y u = \nabla^2\psi$. The equation solved by the module is: -$$\partial_t q + \J(\psi, q) = \underbrace{-\left[\mu(-1)^{n_\mu} \nabla^{2n_\mu} -+\nu(-1)^{n_\nu} \nabla^{2n_\nu}\right] q}_{\textrm{dissipation}} + f\ .$$ +$$\partial_t \zeta + \J(\psi, \zeta) = \underbrace{-\left[\mu(-1)^{n_\mu} \nabla^{2n_\mu} ++\nu(-1)^{n_\nu} \nabla^{2n_\nu}\right] \zeta}_{\textrm{dissipation}} + f\ .$$ where $\J(a, b) = (\partial_x a)(\partial_y b)-(\partial_y a)(\partial_x b)$. On the right hand side, $f(x,y,t)$ is forcing, $\mu$ is hypoviscosity, and $\nu$ is @@ -24,8 +24,8 @@ viscosity corresponds to $n_{\nu}=1$. The equation is time-stepped forward in Fourier space: -$$\partial_t \widehat{q} = - \widehat{J(\psi, q)} -\left(\mu k^{2n_\mu} -+\nu k^{2n_\nu}\right) \widehat{q} + \widehat{f}\ .$$ +$$\partial_t \widehat{\zeta} = - \widehat{J(\psi, \zeta)} -\left(\mu k^{2n_\mu} ++\nu k^{2n_\nu}\right) \widehat{\zeta} + \widehat{f}\ .$$ In doing so the Jacobian is computed in the conservative form: $\J(a,b) = \partial_y [ (\partial_x a) b] -\partial_x[ (\partial_y a) b]$. @@ -33,8 +33,8 @@ In doing so the Jacobian is computed in the conservative form: $\J(a,b) = Thus: $$\mathcal{L} = -\mu k^{-2n_\mu} - \nu k^{2n_\nu}\ ,$$ -$$\mathcal{N}(\widehat{q}) = - \mathrm{i}k_x \mathrm{FFT}(u q)- - \mathrm{i}k_y \mathrm{FFT}(\upsilon q) + \widehat{f}\ .$$ +$$\mathcal{N}(\widehat{\zeta}) = - \mathrm{i}k_x \mathrm{FFT}(u \zeta)- + \mathrm{i}k_y \mathrm{FFT}(\upsilon \zeta) + \widehat{f}\ .$$ ### AbstractTypes and Functions @@ -54,13 +54,13 @@ For the forced case ($f\ne 0$) parameters AbstractType is build with `ForcedPara **Vars** For the unforced case ($f=0$) variables AbstractType is build with `Vars` and it includes: -- `q`: Array of Floats; relative vorticity. -- `U`: Array of Floats; $x$-velocity, $u$. -- `V`: Array of Floats; $y$-velocity, $v$. -- `sol`: Array of Complex; the solution, $\widehat{q}$. -- `qh`: Array of Complex; the Fourier transform $\widehat{q}$. -- `Uh`: Array of Complex; the Fourier transform $\widehat{u}$. -- `Vh`: Array of Complex; the Fourier transform $\widehat{v}$. +- `zeta`: Array of Floats; relative vorticity. +- `u`: Array of Floats; $x$-velocity, $u$. +- `v`: Array of Floats; $y$-velocity, $\upsilon$. +- `sol`: Array of Complex; the solution, $\widehat{\zeta}$. +- `zetah`: Array of Complex; the Fourier transform $\widehat{\zeta}$. +- `uh`: Array of Complex; the Fourier transform $\widehat{u}$. +- `vh`: Array of Complex; the Fourier transform $\widehat{\upsilon}$. For the forced case ($f\ne 0$) variables AbstractType is build with `ForcedVars`. It includes all variables in `Vars` and additionally: - `Fh`: Array of Complex; the Fourier transform $\widehat{f}$. @@ -70,51 +70,21 @@ For the forced case ($f\ne 0$) variables AbstractType is build with `ForcedVars` **`calcN!` function** -The nonlinear term $\mathcal{N}(\widehat{q})$ is computed via functions: +The nonlinear term $\mathcal{N}(\widehat{\zeta})$ is computed via functions: -- `calcN_advection!`: computes $- \widehat{J(\psi, q)}$ and stores it in array `N`. +- `calcN_advection!`: computes $- \widehat{J(\psi, \zeta)}$ and stores it in array `N`. -```julia -function calcN_advection!(N, sol, t, s, v, p, g) - @. v.Uh = im * g.l * g.invKrsq * sol - @. v.Vh = -im * g.kr * g.invKrsq * sol - @. v.qh = sol +- `calcN_forced!`: computes $- \widehat{J(\psi, \zeta)}$ via `calcN_advection!` and then adds to it the forcing $\widehat{f}$ computed via `calcF!` function. Also saves the solution $\widehat{\zeta}$ of the previous time-step in array `prevsol`. - A_mul_B!(v.U, g.irfftplan, v.Uh) - A_mul_B!s(v.V, g.irfftplan, v.Vh) - A_mul_B!(v.q, g.irfftplan, v.qh) - - @. v.U *= v.q # U*q - @. v.V *= v.q # V*q - - A_mul_B!(v.Uh, g.rfftplan, v.U) # \hat{U*q} - A_mul_B!(v.Vh, g.rfftplan, v.V) # \hat{U*q} - - @. N = -im*g.kr*v.Uh - im*g.l*v.Vh - nothing -end -``` - -- `calcN_forced!`: computes $- \widehat{J(\psi, q)}$ via `calcN_advection!` and then adds to it the forcing $\widehat{f}$ computed via `calcF!` function. Also saves the solution $\widehat{q}$ of the previous time-step in array `prevsol`. - -```julia -function calcN_forced!(N, sol, t, s, v, p, g) - calcN_advection!(N, sol, t, s, v, p, g) - if t == s.t # not a substep - v.prevsol .= s.sol # used to compute budgets when forcing is stochastic - p.calcF!(v.Fh, sol, t, s, v, p, g) - end - @. N += v.Fh - nothing -end -``` -- `updatevars!`: uses `sol` to compute $q$, $u$, $v$, $\widehat{u}$, and $\widehat{v}$ and stores them into corresponding arrays of `Vars`/`ForcedVars`. +- `updatevars!`: uses `sol` to compute $\zeta$, $u$, $\upsilon$, $\widehat{u}$, and $\widehat{\upsilon}$ and stores them into corresponding arrays of `Vars`/`ForcedVars`. ## Examples -- `examples/twodturb/McWilliams.jl`: A script that simulates decaying two-dimensional turbulence reproducing the results of the paper by +- `examples/twodturb_mcwilliams1984.jl`: A script that simulates decaying two-dimensional turbulence reproducing the results of the paper by > McWilliams, J. C. (1984). The emergence of isolated coherent vortices in turbulent flow. *J. Fluid Mech.*, **146**, 21-43. -- `examples/twodturb/IsotropicRingForcing.jl`: A script that simulates stochastically forced two-dimensional turbulence. The forcing is temporally delta-corraleted and its spatial structure is isotropic with power in a narrow annulus of total radius $k_f$ in wavenumber space. +- `examples/twodturb_randomdecay.jl`: A script that simulates decaying two-dimensional turbulence starting from random initial conditions. + +- `examples/twodturb_stochasticforcing.jl`: A script that simulates forced-dissipative two-dimensional turbulence with isotropic temporally delta-correlated stochastic forcing. diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 8908971c..89b8b240 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -45,7 +45,7 @@ sol, cl, pr, vs, gr = prob.sol, prob.clock, prob.params, prob.vars, prob.grid filepath = "." plotpath = "./plots_2layer" plotname = "snapshots" -filename = joinpath(filepath, "2layer.jl.jld2") +filename = joinpath(filepath, "2layer.jld2") # File management if isfile(filename); rm(filename); end diff --git a/examples/twodturb_mcwilliams1984.jl b/examples/twodturb_mcwilliams1984.jl index 93f53ee5..a5dc17c3 100644 --- a/examples/twodturb_mcwilliams1984.jl +++ b/examples/twodturb_mcwilliams1984.jl @@ -35,8 +35,8 @@ x, y = gridpoints(gr) # that reproduces the results of the paper by McWilliams (1984) seed!(1234) k0, E0 = 6, 0.5 -qi = peakedisotropicspectrum(gr, k0, E0, mask=filter) -TwoDTurb.set_q!(prob, qi) +zetai = peakedisotropicspectrum(gr, k0, E0, mask=filter) +TwoDTurb.set_zeta!(prob, zetai) # Create Diagnostic -- energy and enstrophy are functions imported at the top. E = Diagnostic(energy, prob; nsteps=nsteps) @@ -45,7 +45,7 @@ diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will # be updated every timestep. # Create Output -get_sol(prob) = prob.vars.sol # extracts the Fourier-transformed solution +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution get_u(prob) = irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx) out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) @@ -54,7 +54,7 @@ function plot_output(prob, fig, axs; drawcolorbar=false) # Plot the vorticity field and the evolution of energy and enstrophy. TwoDTurb.updatevars!(prob) sca(axs[1]) - pcolormesh(x, y, vs.q) + pcolormesh(x, y, vs.zeta) clim(-40, 40) axis("off") axis("square") diff --git a/examples/twodturb_randomdecay.jl b/examples/twodturb_randomdecay.jl index a715f5e3..8a82db05 100644 --- a/examples/twodturb_randomdecay.jl +++ b/examples/twodturb_randomdecay.jl @@ -18,7 +18,7 @@ ntot = 10nint # Number of total timesteps # Define problem prob = TwoDTurb.Problem(nx=nx, Lx=Lx, nu=nu, nnu=nnu, dt=dt, stepper="FilteredRK4") -TwoDTurb.set_q!(prob, rand(nx, nx)) +TwoDTurb.set_zeta!(prob, rand(nx, nx)) cl, vs, gr = prob.clock, prob.vars, prob.grid x, y = gridpoints(gr) @@ -27,7 +27,7 @@ x, y = gridpoints(gr) function makeplot!(ax, prob) sca(ax) cla() - pcolormesh(x, y, vs.q) + pcolormesh(x, y, vs.zeta) title("Vorticity") xlabel(L"x") ylabel(L"y") @@ -56,7 +56,7 @@ kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1) fig2, axs = subplots(ncols=2, figsize=(8, 4)) sca(axs[1]) -pcolormesh(x, y, vs.q) +pcolormesh(x, y, vs.zeta) xlabel(L"x") ylabel(L"y") title("Vorticity") diff --git a/examples/twodturb_stochasticforcing.jl b/examples/twodturb_stochasticforcing.jl index 47eb5cea..687f7bdb 100644 --- a/examples/twodturb_stochasticforcing.jl +++ b/examples/twodturb_stochasticforcing.jl @@ -43,7 +43,7 @@ prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, nmu=nmu, dt=dt, stepp sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -TwoDTurb.set_q!(prob, 0*x) +TwoDTurb.set_zeta!(prob, 0*x) E = Diagnostic(energy, prob, nsteps=nt) D = Diagnostic(dissipation, prob, nsteps=nt) R = Diagnostic(drag, prob, nsteps=nt) @@ -57,7 +57,7 @@ function makeplot(prob, diags) t = round(mu*cl.t, digits=2) sca(axs[1]); cla() - pcolormesh(x, y, v.q) + pcolormesh(x, y, v.zeta) xlabel(L"$x$") ylabel(L"$y$") title("\$\\nabla^2\\psi(x,y,\\mu t= $t )\$") diff --git a/src/twodturb.jl b/src/twodturb.jl index f5f8d732..7b6e9451 100644 --- a/src/twodturb.jl +++ b/src/twodturb.jl @@ -2,7 +2,7 @@ module TwoDTurb export Problem, - set_q!, + set_zeta!, updatevars!, energy, @@ -22,7 +22,7 @@ using FourierFlows: getfieldspecs, structvarsexpr, parsevalsum, parsevalsum2 abstract type TwoDTurbVars <: AbstractVars end -const physicalvars = [:q, :u, :v] +const physicalvars = [:zeta, :u, :v] const transformvars = [ Symbol(var, :h) for var in physicalvars ] const forcedvars = [:Fh] const stochforcedvars = [:prevsol] @@ -90,9 +90,9 @@ Params(nu, nnu) = Params(nu, nnu, typeof(nu)(0), 0, nothingfunction) Returns the equation for two-dimensional turbulence with params p and grid g. """ function Equation(p::Params, g::AbstractGrid{T}) where T - LC = @. - p.nu*g.Krsq^p.nnu - p.mu*g.Krsq^p.nmu - LC[1, 1] = 0 - FourierFlows.Equation(LC, calcN!, g) + L = @. - p.nu*g.Krsq^p.nnu - p.mu*g.Krsq^p.nmu + L[1, 1] = 0 + FourierFlows.Equation(L, calcN!, g) end @@ -119,9 +119,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) zeta u v + @createarrays Complex{T} (g.nkr, g.nl) sol zetah uh vh + Vars(zeta, u, v, zetah, uh, vh) end """ @@ -160,17 +160,17 @@ 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.qh = sol + @. v.zetah = sol ldiv!(v.u, g.rfftplan, v.uh) ldiv!(v.v, g.rfftplan, v.vh) - ldiv!(v.q, g.rfftplan, v.qh) + ldiv!(v.zeta, g.rfftplan, v.zetah) - @. v.u *= v.q # U*q - @. v.v *= v.q # V*q + @. v.u *= v.zeta # u*zeta + @. v.v *= v.zeta # v*zeta - 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*zeta} + mul!(v.vh, g.rfftplan, v.v) # \hat{v*zeta} @. N = -im*g.kr*v.uh - im*g.l*v.vh nothing @@ -211,24 +211,24 @@ 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.zetah .= 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.zeta, g.rfftplan, deepcopy(v.zetah)) ldiv!(v.u, g.rfftplan, deepcopy(v.uh)) ldiv!(v.v, g.rfftplan, deepcopy(v.vh)) nothing end """ - set_q!(prob, q) + set_zeta!(prob, zeta) -Set the solution sol as the transform of q and update variables v +Set the solution sol as the transform of zeta and update variables v on the grid g. """ -function set_q!(prob, q) +function set_zeta!(prob, zeta) p, v, g, sol = prob.params, prob.vars, prob.grid, prob.sol - mul!(sol, g.rfftplan, q) + mul!(sol, g.rfftplan, zeta) sol[1, 1] = 0 # zero domain average updatevars!(prob) nothing diff --git a/test/test_twodturb.jl b/test/test_twodturb.jl index 3b61a77d..ec73a09b 100644 --- a/test/test_twodturb.jl +++ b/test/test_twodturb.jl @@ -1,20 +1,20 @@ function test_twodturb_lambdipole(n, dt; L=2π, Ue=1, Re=L/20, nu=0.0, nnu=1, ti=L/Ue*0.01, nm=3) nt = round(Int, ti/dt) prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, dt=dt, stepper="FilteredRK4") - q₀ = lambdipole(Ue, Re, prob.grid) - TwoDTurb.set_q!(prob, q₀) + zeta₀ = lambdipole(Ue, Re, prob.grid) + TwoDTurb.set_zeta!(prob, zeta₀) - xq = zeros(nm) # centroid of abs(q) + xzeta = zeros(nm) # centroid of abs(zeta) Ue_m = zeros(nm) # measured dipole speed x, y = gridpoints(prob.grid) - q = prob.vars.q + zeta = prob.vars.zeta for i = 1:nm # step forward stepforward!(prob, nt) TwoDTurb.updatevars!(prob) - xq[i] = mean(abs.(q).*x) / mean(abs.(q)) + xzeta[i] = mean(@. abs(zeta)*x) / mean(abs.(zeta)) if i > 1 - Ue_m[i] = (xq[i]-xq[i-1]) / ((nt-1)*dt) + Ue_m[i] = (xzeta[i]-xzeta[i-1]) / ((nt-1)*dt) end end isapprox(Ue, mean(Ue_m[2:end]), rtol=rtol_lambdipole) @@ -54,7 +54,7 @@ function test_twodturb_stochasticforcingbudgets(; n=256, L=2π, dt=0.005, nu=1e- sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid; - TwoDTurb.set_q!(prob, 0*x) + TwoDTurb.set_zeta!(prob, 0*x) E = Diagnostic(TwoDTurb.energy, prob, nsteps=nt) D = Diagnostic(TwoDTurb.dissipation, prob, nsteps=nt) R = Diagnostic(TwoDTurb.drag, prob, nsteps=nt) @@ -108,7 +108,7 @@ function test_twodturb_deterministicforcingbudgets(; n=256, dt=0.01, L=2π, nu=1 sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid - TwoDTurb.set_q!(prob, 0*x) + TwoDTurb.set_zeta!(prob, 0*x) E = Diagnostic(TwoDTurb.energy, prob, nsteps=nt) D = Diagnostic(TwoDTurb.dissipation, prob, nsteps=nt) @@ -156,8 +156,8 @@ function test_twodturb_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0 gr = TwoDGrid(n, L) x, y = gridpoints(gr) - psif = @. sin(2x)*cos(2y) + 2sin(x)*cos(3y) - qf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y) + psif = @. sin(2x)*cos(2y) + 2sin(x)*cos(3y) + zetaf = @. -8sin(2x)*cos(2y) - 20sin(x)*cos(3y) Ff = @. -( nu*( 64sin(2x)*cos(2y) + 200sin(x)*cos(3y) ) @@ -174,12 +174,12 @@ function test_twodturb_advection(dt, stepper; n=128, L=2π, nu=1e-2, nnu=1, mu=0 prob = TwoDTurb.Problem(nx=n, Lx=L, nu=nu, nnu=nnu, mu=mu, nmu=nmu, dt=dt, stepper=stepper, calcF=calcF!, stochastic=false) sol, cl, p, v, g = prob.sol, prob.clock, prob.params, prob.vars, prob.grid - TwoDTurb.set_q!(prob, qf) + TwoDTurb.set_zeta!(prob, zetaf) stepforward!(prob, nt) TwoDTurb.updatevars!(prob) - isapprox(prob.vars.q, qf, rtol=rtol_twodturb) + isapprox(prob.vars.zeta, zetaf, rtol=rtol_twodturb) end function test_twodturb_energyenstrophy() @@ -189,20 +189,20 @@ function test_twodturb_energyenstrophy() k0, l0 = gr.k[2], gr.l[2] # fundamental wavenumbers x, y = gridpoints(gr) - psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y) - q0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y) + psi0 = @. sin(2k0*x)*cos(2l0*y) + 2sin(k0*x)*cos(3l0*y) + zeta0 = @. -((2k0)^2+(2l0)^2)*sin(2k0*x)*cos(2l0*y) - (k0^2+(3l0)^2)*2sin(k0*x)*cos(3l0*y) energy_calc = 29/9 enstrophy_calc = 2701/162 prob = TwoDTurb.Problem(nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler") - TwoDTurb.set_q!(prob, q0) + TwoDTurb.set_zeta!(prob, zeta0) TwoDTurb.updatevars!(prob) - energyq0 = TwoDTurb.energy(prob) - enstrophyq0 = TwoDTurb.enstrophy(prob) + energyzeta0 = TwoDTurb.energy(prob) + enstrophyzeta0 = TwoDTurb.enstrophy(prob) - (isapprox(energyq0, 29.0/9, rtol=rtol_twodturb) && - isapprox(enstrophyq0, 2701.0/162, rtol=rtol_twodturb)) + (isapprox(energyzeta0, 29.0/9, rtol=rtol_twodturb) && + isapprox(enstrophyzeta0, 2701.0/162, rtol=rtol_twodturb)) end