Merge bdb85b2 into 18feb28

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,17 +24,17 @@ 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]$.
 
 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.

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

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")

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")

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 )\$")

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

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