Merge pull request #108 from FourierFlows/tweak2DNSstochasticforcinge…

…xample

Separate 2D Navier-Stokes stochastic forcing & budgets examples

docs/make.jl

@@ -21,6 +21,7 @@ const OUTPUT_DIR   = joinpath(@__DIR__, "src/generated")
 examples = [
     "twodnavierstokes_decaying.jl",
     "twodnavierstokes_stochasticforcing.jl",
+    "twodnavierstokes_stochasticforcing_budgets.jl",
     "barotropicqg_betadecay.jl",
     "barotropicqg_betaforced.jl",
     "barotropicqg_acc.jl",
@@ -73,6 +74,7 @@ sitename = "GeophysicalFlows.jl",
               "TwoDNavierStokes" => Any[
                 "generated/twodnavierstokes_decaying.md",
                 "generated/twodnavierstokes_stochasticforcing.md",
+                "generated/twodnavierstokes_stochasticforcing_budgets.md",
                 ],
               "BarotropicQG" => Any[
                 "generated/barotropicqg_betadecay.md",

examples/barotropicqg_betadecay.jl

@@ -22,7 +22,7 @@ nothing # hide
 
 # ## Numerical parameters and time-stepping parameters
 
-     nx = 128            # 2D resolution = nx^2
+      n = 128            # 2D resolution = n^2
 stepper = "FilteredRK4"  # timestepper
      dt = 0.05           # timestep
  nsteps = 2000           # total number of time-steps
@@ -32,17 +32,17 @@ nothing # hide
 
 # ## Physical parameters
 
-Lx = 2π        # domain size
- β = 10.0      # planetary PV gradient
- μ = 0.0       # bottom drag
+L = 2π        # domain size
+β = 10.0      # planetary PV gradient
+μ = 0.0       # bottom drag
 nothing # hide
 
 # ## Problem setup
 # We initialize a `Problem` by providing a set of keyword arguments. Not providing
 # a viscosity coefficient ν leads to the module's default value: ν=0. In this
 # example numerical instability due to accumulation of enstrophy in high wavenumbers
 # is taken care with the `FilteredTimestepper` we picked. 
-prob = BarotropicQG.Problem(dev; nx=nx, Lx=Lx, β=β, μ=μ, dt=dt, stepper=stepper)
+prob = BarotropicQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper)
 nothing # hide
 
 # and define some shortcuts
@@ -54,19 +54,20 @@ nothing # hide
 # ## Setting initial conditions
 
 # Our initial condition consist of a flow that has power only at wavenumbers with
-# $8<\frac{L}{2\pi}\sqrt{k_x^2+k_y^2}<10$ and initial energy $E_0$:
+# $8 < \frac{L}{2\pi} \sqrt{k_x^2+k_y^2} < 10$ and initial energy $E_0$:
 
 E0 = 0.1 # energy of initial condition
 
-k = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl] # a 2D grid with the zonal wavenumber
+K = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
+k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]  # a 2D array with the zonal wavenumber
 
 Random.seed!(1234)
 qih = randn(Complex{eltype(gr)}, size(sol))
-qih[ gr.Krsq .< (8*2π /gr.Lx)^2 ] .= 0
-qih[ gr.Krsq .> (10*2π/gr.Lx)^2 ] .= 0
-qih[ k .== 0 ] .= 0 # remove any power from k=0 component
-Ein = energy(qih, gr)  # compute energy of qi
-qih = qih*sqrt(E0/Ein) # normalize qi to have energy E0
+qih[K .<  8 * 2π/L] .= 0
+qih[K .> 10 * 2π/L] .= 0
+qih[k .== 0] .= 0         # no power at zonal wavenumber k=0 component
+Ein = energy(qih, gr)     # compute energy of qi
+qih *= sqrt(E0 / Ein)     # normalize qi to have energy E0
 qi  = irfft(qih, gr.nx)
 
 BarotropicQG.set_zeta!(prob, qi)

examples/barotropicqg_betaforced.jl

@@ -25,7 +25,7 @@ nothing # hide
 
 # ## Numerical parameters and time-stepping parameters
 
-     nx = 128            # 2D resolution = nx^2
+      n = 128            # 2D resolution = n^2
 stepper = "FilteredRK4"  # timestepper
      dt = 0.05           # timestep
  nsteps = 8000           # total number of time-steps
@@ -35,9 +35,9 @@ nothing # hide
 
 # ## Physical parameters
 
-Lx = 2π        # domain size
- β = 10.0      # planetary PV gradient
- μ = 0.01      # bottom drag
+L = 2π        # domain size
+β = 10.0      # planetary PV gradient
+μ = 0.01      # bottom drag
 nothing # hide
 
 
@@ -53,26 +53,27 @@ forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum tha
 forcing_bandwidth  = 1.5     # the width of the forcing spectrum 
 ε = 0.001                    # energy input rate by the forcing
 
-gr  = TwoDGrid(nx, Lx)
+gr  = TwoDGrid(n, L)
 
-k = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl] # a 2D grid with the zonal wavenumber
+K  = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
+k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]   # a 2D array with the zonal wavenumber
 
-forcing_spectrum = @. exp( -(sqrt(gr.Krsq)-forcing_wavenumber)^2 / (2forcing_bandwidth^2) )
-@. forcing_spectrum[ gr.Krsq < (2π/Lx*2)^2  ] = 0
-@. forcing_spectrum[ gr.Krsq > (2π/Lx*20)^2 ] = 0
-@. forcing_spectrum[ k .< 2π/Lx ] .= 0 # make sure forcing does not have power at k=0 component
-ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq/2, gr)/(gr.Lx*gr.Ly)
-@. forcing_spectrum = ε/ε0 * forcing_spectrum  # normalization so that forcing injects energy ε per domain area per unit time
+forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
+forcing_spectrum[K .<  2 * 2π/L] .= 0    # no power at low wavenumbers
+forcing_spectrum[K .> 20 * 2π/L] .= 0    # no power at high wavenumbers
+forcing_spectrum[k .< 2π/L] .= 0         # make sure forcing does not have power at k=0
+ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
+@. forcing_spectrum *= ε/ε0               # normalize forcing to inject energy at rate ε
 
 seed!(1234) # reset of the random number generator for reproducibility
 nothing # hide
 
 # Next we construct function `calcF!` that computes a forcing realization every timestep
 function calcFq!(Fh, sol, t, clock, vars, params, grid)
-  ξ = ArrayType(dev)(exp.(2π*im*rand(eltype(grid), size(sol)))/sqrt(clock.dt))
+  ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
   @. Fh = ξ*sqrt.(forcing_spectrum)
-  Fh[abs.(grid.Krsq).==0] .= 0
-  nothing
+  Fh[abs.(grid.Krsq) .== 0] .= 0
+  return nothing
 end
 nothing # hide
 
@@ -82,7 +83,7 @@ nothing # hide
 # a viscosity coefficient ν leads to the module's default value: ν=0. In this
 # example numerical instability due to accumulation of enstrophy in high wavenumbers
 # is taken care with the `FilteredTimestepper` we picked. 
-prob = BarotropicQG.Problem(dev; nx=nx, Lx=Lx, β=β, μ=μ, dt=dt, stepper=stepper, 
+prob = BarotropicQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper, 
                             calcFq=calcFq!, stochastic=true)
 nothing # hide
 

examples/barotropicqgql_betaforced.jl

@@ -26,7 +26,7 @@ nothing # hide
 
 # ## Numerical parameters and time-stepping parameters
 
-     nx = 128            # 2D resolution = nx^2
+      n = 128            # 2D resolution = n^2
 stepper = "FilteredRK4"  # timestepper
      dt = 0.05           # timestep
  nsteps = 8000           # total number of time-steps
@@ -36,9 +36,9 @@ nothing # hide
 
 # ## Physical parameters
 
-Lx = 2π        # domain size
- β = 10.0      # planetary PV gradient
- μ = 0.01      # bottom drag
+L = 2π        # domain size
+β = 10.0      # planetary PV gradient
+μ = 0.01      # bottom drag
 nothing # hide
 
 
@@ -54,26 +54,27 @@ forcing_wavenumber = 14.0    # the central forcing wavenumber for a spectrum tha
 forcing_bandwidth  = 1.5     # the width of the forcing spectrum 
 ε = 0.001                    # energy input rate by the forcing
 
-gr  = TwoDGrid(nx, Lx)
+gr  = TwoDGrid(n, L)
 
-k = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl ] # a 2D grid with the zonal wavenumber
+K = @. sqrt(gr.Krsq)                      # a 2D array with the total wavenumber
+k = [gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]  # a 2D array with the zonal wavenumber
 
-forcing_spectrum = @. exp( -(sqrt(gr.Krsq)-forcing_wavenumber)^2 / (2forcing_bandwidth^2) )
-@. forcing_spectrum[ gr.Krsq < (2π/Lx*2)^2  ] = 0
-@. forcing_spectrum[ gr.Krsq > (2π/Lx*20)^2 ] = 0
-@. forcing_spectrum[ k .< 2π/Lx ] .= 0 # make sure forcing does not have power at k=0 component
-ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq/2, gr) / (gr.Lx*gr.Ly)
-@. forcing_spectrum = ε/ε0 * forcing_spectrum  # normalization so that forcing injects energy ε per domain area per unit time
+forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2))
+forcing_spectrum[K .<  2 * 2π/L] .= 0    # no power at low wavenumbers
+forcing_spectrum[K .> 20 * 2π/L] .= 0    # no power at high wavenumbers
+forcing_spectrum[k .< 2π/L ] .= 0        # make sure forcing does not have power at k=0
+ε0 = parsevalsum(forcing_spectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly)
+@. forcing_spectrum *= ε/ε0               # normalize forcing to inject energy at rate ε
 
 seed!(1234) # reset of the random number generator for reproducibility
 nothing # hide
 
 # Next we construct function `calcF!` that computes a forcing realization every timestep
 function calcF!(Fh, sol, t, clock, vars, params, grid)
-  ξ = ArrayType(dev)(exp.(2π*im*rand(eltype(grid), size(sol)))/sqrt(clock.dt))
-  @. Fh = ξ*sqrt.(forcing_spectrum)
-  Fh[abs.(grid.Krsq).==0] .= 0
-  nothing
+  ξ = ArrayType(dev)(exp.(2π * im * rand(eltype(grid), size(sol))) / sqrt(clock.dt))
+  @. Fh = ξ * sqrt.(forcing_spectrum)
+  Fh[abs.(grid.Krsq) .== 0] .= 0
+  return nothing
 end
 nothing # hide
 
@@ -83,7 +84,7 @@ nothing # hide
 # a viscosity coefficient ν leads to the module's default value: ν=0. In this
 # example numerical instability due to accumulation of enstrophy in high wavenumbers
 # is taken care with the `FilteredTimestepper` we picked. 
-prob = BarotropicQGQL.Problem(dev; nx=nx, Lx=Lx, β=β, μ=μ, dt=dt, stepper=stepper, 
+prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper, 
                               calcF=calcF!, stochastic=true)
 nothing # hide
 

examples/multilayerqg_2layer.jl

@@ -26,9 +26,9 @@ nx = 128        # 2D resolution = nx^2
 ny = nx
 
 stepper = "FilteredRK4"  # timestepper
-dt = 6e-3                # timestep
-nsteps = 7000            # total number of time-steps
-nsubs  = 25              # number of time-steps for plotting (nsteps must be multiple of nsubs)
+     dt = 6e-3           # timestep
+ nsteps = 7000           # total number of time-steps
+ nsubs  = 25             # number of time-steps for plotting (nsteps must be multiple of nsubs)
 nothing # hide