From 43411a6380299378166ac32db217da5542924ccc Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 08:52:38 +1000 Subject: [PATCH 1/9] clears up the 2DNS budgets example --- docs/make.jl | 4 +- .../twodnavierstokes_stochasticforcing.jl | 248 ------------------ 2 files changed, 2 insertions(+), 250 deletions(-) delete mode 100644 examples/twodnavierstokes_stochasticforcing.jl diff --git a/docs/make.jl b/docs/make.jl index 49b28f3f..a96d1f4c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,7 +20,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", @@ -72,7 +72,7 @@ sitename = "GeophysicalFlows.jl", "Examples" => [ "TwoDNavierStokes" => Any[ "generated/twodnavierstokes_decaying.md", - "generated/twodnavierstokes_stochasticforcing.md", + "generated/twodnavierstokes_stochasticforcing_budgets.md", ], "BarotropicQG" => Any[ "generated/barotropicqg_betadecay.md", diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl deleted file mode 100644 index 9468c47b..00000000 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ /dev/null @@ -1,248 +0,0 @@ -# # 2D forced-dissipative turbulence -# -#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). -#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). -# -# A simulation of forced-dissipative two-dimensional turbulence. We solve the -# two-dimensional vorticity equation with stochastic excitation and dissipation in -# the form of linear drag and hyperviscosity. As a demonstration, we compute how -# each of the forcing and dissipation terms contribute to the energy and the -# enstrophy budgets. - -using FourierFlows, Printf, Plots - -using FourierFlows: parsevalsum -using Random: seed! -using FFTW: irfft - -import GeophysicalFlows.TwoDNavierStokes -import GeophysicalFlows.TwoDNavierStokes: energy, energy_dissipation, energy_work, energy_drag -import GeophysicalFlows.TwoDNavierStokes: enstrophy, enstrophy_dissipation, enstrophy_work, enstrophy_drag - - -# ## Choosing a device: CPU or GPU - -dev = CPU() # Device (CPU/GPU) -nothing # hide - - -# ## Numerical, domain, and simulation parameters -# -# First, we pick some numerical and physical parameters for our model. - - n, L = 256, 2π # grid resolution and domain length - ν, nν = 2e-7, 2 # hyperviscosity coefficient and order - μ, nμ = 1e-1, 0 # linear drag coefficient -dt, tf = 0.005, 0.2/μ # timestep and final time - nt = round(Int, tf/dt) # total timesteps - ns = 4 # how many intermediate times we want to plot -nothing # hide - - -# ## Forcing - -# We force the vorticity equation with stochastic excitation that is delta-correlated -# in time and while spatially homogeneously and isotropically correlated. The forcing -# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and -# width $\delta k_f$, and it injects energy per unit area and per unit time equal -# to $\varepsilon$. - -forcing_wavenumber = 14.0 # the central forcing wavenumber for a spectrum that is a ring in wavenumber space -forcing_bandwidth = 1.5 # the width of the forcing spectrum -ε = 0.001 # energy input rate by the forcing - -gr = TwoDGrid(dev, n, L) -x, y = gr.x, gr.y - -forcing_spectrum = @. exp(-(sqrt(gr.Krsq)-forcing_wavenumber)^2/(2*forcing_bandwidth^2)) -forcing_spectrum[ gr.Krsq .< (2π/L*2)^2 ] .= 0 # make sure that focing has no power for low wavenumbers -forcing_spectrum[ gr.Krsq .> (2π/L*20)^2 ] .= 0 # make sure that focing has no power for high wavenumbers -ε0 = parsevalsum(forcing_spectrum.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly) -forcing_spectrum .= ε/ε0 * forcing_spectrum # normalize forcing to inject energy ε - -seed!(1234) -nothing # hide - -# Next we construct function `calcF!` that computes a forcing realization every timestep -function calcF!(Fh, sol, t, cl, v, p, g) - ξ = ArrayType(dev)(exp.(2π*im*rand(eltype(gr), size(sol)))/sqrt(cl.dt)) - ξ[1, 1] = 0 - @. Fh = ξ*sqrt(forcing_spectrum) - nothing -end -nothing # hide - - -# ## Problem setup -# We initialize a `Problem` by providing a set of keyword arguments. The -# `stepper` keyword defines the time-stepper to be used. -prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4", - calcF=calcF!, stochastic=true) -nothing # hide - -# Define some shortcuts for convenience. -sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -nothing # hide - - -# First let's see how a forcing realization looks like. -calcF!(v.Fh, sol, 0.0, cl, v, p, g) - -heatmap(x, y, irfft(v.Fh, g.nx), - aspectratio = 1, - c = :balance, - clim = (-200, 200), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "a forcing realization", - framestyle = :box) - - -# ## Setting initial conditions - -# Our initial condition is simply fluid at rest. -TwoDNavierStokes.set_zeta!(prob, zeros(g.nx, g.ny)) - - -# ## Diagnostics - -# Create Diagnostics; the diagnostics are aimed to probe the energy budget. -E = Diagnostic(energy, prob, nsteps=nt) # energy -R = Diagnostic(energy_drag, prob, nsteps=nt) # dissipation by drag -D = Diagnostic(energy_dissipation, prob, nsteps=nt) # dissipation by hyperviscosity -W = Diagnostic(energy_work, prob, nsteps=nt) # work input by forcing -Z = Diagnostic(enstrophy, prob, nsteps=nt) # energy -RZ = Diagnostic(enstrophy_drag, prob, nsteps=nt) # dissipation by drag -DZ = Diagnostic(enstrophy_dissipation, prob, nsteps=nt) # dissipation by hyperviscosity -WZ = Diagnostic(enstrophy_work, prob, nsteps=nt) # work input by forcing -diags = [E, D, W, R, Z, DZ, WZ, RZ] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. -nothing # hide - - -# ## Visualizing the simulation - -# We define a function that plots the vorticity field and the evolution of -# the diagnostics: energy and all terms involved in the energy budget. Last -# we confirm whether the energy budget is accurate, i.e., $\mathrm{d}E/\mathrm{d}t = W - R - D$. - -function computetendencies_and_makeplot(prob, diags) - TwoDNavierStokes.updatevars!(prob) - E, D, W, R, Z, Dᶻ, Wᶻ, Rᶻ = diags - - clocktime = round(μ*cl.t, digits=2) - - i₀ = 1 - dEdt_numerical = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt #numerical first-order approximation of energy tendency - dZdt_numerical = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt #numerical first-order approximation of enstrophy tendency - ii = (i₀):E.i-1 - ii2 = (i₀+1):E.i - - t = E.t[ii] - dEdt_computed = W[ii2] - D[ii] - R[ii] # Stratonovich interpretation - dZdt_computed = Wᶻ[ii2] - Dᶻ[ii] - Rᶻ[ii] - - residual_E = dEdt_computed - dEdt_numerical - residual_Z = dZdt_computed - dZdt_numerical - - εᶻ = parsevalsum(forcing_spectrum / 2, g) / (g.Lx * g.Ly) - - l = @layout grid(2,3) - - pzeta = heatmap(x, y, v.zeta, - aspectratio = 1, - legend = false, - c = :viridis, - clim = (-25, 25), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "μt", - ylabel = "y", - title = "∇²ψ(x, y, t="*@sprintf("%.2f", cl.t)*")", - framestyle = :box) - - pζ = plot(pzeta, size = (400, 400)) - - p1 = plot(μ*t, [W[ii2] ε.+0*t -D[ii] -R[ii]], - label = ["work, W" "ensemble mean work, " "dissipation, D" "drag, R=-2μE"], - linestyle = [:solid :dash :solid :solid], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "energy sources and sinks") - - p2 = plot(μ*t, [dEdt_computed[ii], dEdt_numerical], - label = ["computed W-D" "numerical dE/dt"], - linestyle = [:solid :dashdotdot], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "dE/dt") - - p3 = plot(μ*t, residual_E, - label = "residual dE/dt = computed - numerical", - linewidth = 2, - alpha = 0.7, - xlabel = "μt") - - p4 = plot(μ*t, [Wᶻ[ii2] εᶻ.+0*t -Dᶻ[ii] -Rᶻ[ii]], - label = ["Enstrophy work, Wᶻ" "mean enstrophy work, " "enstrophy dissipation, Dᶻ" "enstrophy drag, Rᶻ=-2μZ"], - linestyle = [:solid :dash :solid :solid], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "enstrophy sources and sinks") - - - p5 = plot(μ*t, [dZdt_computed[ii], dZdt_numerical], - label = ["computed Wᶻ-Dᶻ" "numerical dZ/dt"], - linestyle = [:solid :dashdotdot], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "dZ/dt") - - p6 = plot(μ*t, residual_Z, - label = "residual dZ/dt = computed - numerical", - linewidth = 2, - alpha = 0.7, - xlabel = "μt") - - - pbudgets = plot(p1, p2, p3, p4, p5, p6, layout=l, size = (1300, 900)) - - return pζ, pbudgets -end -nothing # hide - - - - -# ## Time-stepping the `Problem` forward - -# Finally, we time-step the `Problem` forward in time. - -startwalltime = time() -for i = 1:ns - stepforward!(prob, diags, round(Int, nt/ns)) - TwoDNavierStokes.updatevars!(prob) - cfl = cl.dt*maximum([maximum(v.u)/g.dx, maximum(v.v)/g.dy]) - - log = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min", cl.step, cl.t, - cfl, (time()-startwalltime)/60) - - println(log) -end - - -# ## Plot -# And now let's see what we got. We plot the output. - -pζ, pbudgets = computetendencies_and_makeplot(prob, diags) - -pbudgets From 51d3779d22429ec2712c1490ac046c03557470b8 Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 08:52:51 +1000 Subject: [PATCH 2/9] clears up the 2DNS budgets example --- ...dnavierstokes_stochasticforcing_budgets.jl | 260 ++++++++++++++++++ 1 file changed, 260 insertions(+) create mode 100644 examples/twodnavierstokes_stochasticforcing_budgets.jl diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl new file mode 100644 index 00000000..ce832f63 --- /dev/null +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -0,0 +1,260 @@ +# # 2D forced-dissipative turbulence +# +#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +# +# A simulation of forced-dissipative two-dimensional turbulence. We solve the +# two-dimensional vorticity equation with stochastic excitation and dissipation in +# the form of linear drag and hyperviscosity. As a demonstration, we compute how +# each of the forcing and dissipation terms contribute to the energy and the +# enstrophy budgets. + +using FourierFlows, Printf, Plots + +using FourierFlows: parsevalsum +using Random: seed! +using FFTW: irfft + +import GeophysicalFlows.TwoDNavierStokes +import GeophysicalFlows.TwoDNavierStokes: energy, energy_dissipation, energy_work, energy_drag +import GeophysicalFlows.TwoDNavierStokes: enstrophy, enstrophy_dissipation, enstrophy_work, enstrophy_drag + + +# ## Choosing a device: CPU or GPU + +dev = CPU() # Device (CPU/GPU) +nothing # hide + + +# ## Numerical, domain, and simulation parameters +# +# First, we pick some numerical and physical parameters for our model. + + n, L = 256, 2π # grid resolution and domain length + ν, nν = 2e-7, 2 # hyperviscosity coefficient and order + μ, nμ = 1e-1, 0 # linear drag coefficient +dt, tf = 0.005, 0.2 / μ # timestep and final time + nt = round(Int, tf / dt) # total timesteps + ns = 4 # how many intermediate times we want to plot +nothing # hide + + +# ## Forcing + +# We force the vorticity equation with stochastic excitation that is delta-correlated +# in time and while spatially homogeneously and isotropically correlated. The forcing +# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and +# width $\delta k_f$, and it injects energy per unit area and per unit time equal +# to $\varepsilon$. + +forcing_wavenumber = 14.0 # the central forcing wavenumber for a spectrum that is a ring in wavenumber space +forcing_bandwidth = 1.5 # the width of the forcing spectrum +ε = 0.1 # energy input rate by the forcing + +grid = TwoDGrid(dev, n, L) +x, y = grid.x, grid.y + +K = sqrt.(grid.Krsq) +forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2)) +forcing_spectrum[K .< (2*2π/L)] .= 0 # make sure that focing has no power for low wavenumbers +forcing_spectrum[K .> (20*2π/L)] .= 0 # make sure that focing has no power for high wavenumbers +ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) +@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy ε + +seed!(1234) +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)) + ξ[1, 1] = 0 + @. Fh = ξ * sqrt(forcing_spectrum) + nothing +end +nothing # hide + + +# ## Problem setup +# We initialize a `Problem` by providing a set of keyword arguments. The +# `stepper` keyword defines the time-stepper to be used. +prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4", + calcF=calcF!, stochastic=true) +nothing # hide + +# Define some shortcuts for convenience. +sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid +nothing # hide + + +# First let's see how a forcing realization looks like. +calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) + +heatmap(x, y, irfft(vars.Fh, grid.nx), + aspectratio = 1, + c = :balance, + clim = (-200, 200), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "x", + ylabel = "y", + title = "a forcing realization", + framestyle = :box) + + +# ## Setting initial conditions + +# Our initial condition is simply fluid at rest. +TwoDNavierStokes.set_zeta!(prob, zeros(grid.nx, grid.ny)) + + +# ## Diagnostics + +# Create Diagnostics; the diagnostics are aimed to probe the energy budget. +E = Diagnostic(energy, prob, nsteps=nt) # energy +Rᵋ = Diagnostic(energy_drag, prob, nsteps=nt) # energy dissipation by drag +Dᵋ = Diagnostic(energy_dissipation, prob, nsteps=nt) # energy dissipation by hyperviscosity +Wᵋ = Diagnostic(energy_work, prob, nsteps=nt) # energy work input by forcing +Z = Diagnostic(enstrophy, prob, nsteps=nt) # enstrophy +Rᶻ = Diagnostic(enstrophy_drag, prob, nsteps=nt) # enstrophy dissipation by drag +Dᶻ = Diagnostic(enstrophy_dissipation, prob, nsteps=nt) # enstrophy dissipation by hyperviscosity +Wᶻ = Diagnostic(enstrophy_work, prob, nsteps=nt) # enstrophy work input by forcing +diags = [E, Dᵋ, Wᵋ, Rᵋ, Z, Dᶻ, Wᶻ, Rᶻ] # a list of Diagnostics passed to `stepforward!` will be updated every timestep. +nothing # hide + + +# ## Visualizing the simulation + +# We define a function that plots the vorticity field and the evolution of +# the diagnostics: energy, enstrophy and all terms involved in the energy and +# enstrophy budgets. Last, we confirm whether the energy and enstrophy budgets +# are accurately computed, e.g., $\mathrm{d}E/\mathrm{d}t = W^\varepsilon - +# R^\varepsilon - D^\varepsilon$. + +function computetendencies_and_makeplot(prob, diags) + sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + + TwoDNavierStokes.updatevars!(prob) + + E, Dᵋ, Wᵋ, Rᵋ, Z, Dᶻ, Wᶻ, Rᶻ = diags + + clocktime = round(μ * clock.t, digits=2) + + dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / clock.dt # numerical first-order approximation of energy tendency + dZdt_numerical = (Z[2:Z.i] - Z[1:Z.i-1]) / clock.dt # numerical first-order approximation of enstrophy tendency + + dEdt_computed = Wᵋ[2:E.i] - Dᵋ[1:E.i-1] - Rᵋ[1:E.i-1] + dZdt_computed = Wᶻ[2:Z.i] - Dᶻ[1:Z.i-1] - Rᶻ[1:Z.i-1] + + residual_E = dEdt_computed - dEdt_numerical + residual_Z = dZdt_computed - dZdt_numerical + + εᶻ = parsevalsum(forcing_spectrum / 2, grid) / (grid.Lx * grid.Ly) + + pzeta = heatmap(x, y, vars.zeta, + aspectratio = 1, + legend = false, + c = :viridis, + clim = (-25, 25), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "μt", + ylabel = "y", + title = "∇²ψ(x, y, μt="*@sprintf("%.2f", μ*clock.t)*")", + framestyle = :box) + + pζ = plot(pzeta, size = (400, 400)) + + t = E.t[2:E.i] + + p1E = plot(μ * t, [Wᵋ[2:E.i] ε.+0*t -Dᵋ[1:E.i-1] -Rᵋ[1:E.i-1]], + label = ["energy work, Wᵋ" "ensemble mean energy work, " "dissipation, Dᵋ" "drag, Rᵋ = - 2μE"], + linestyle = [:solid :dash :solid :solid], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "energy sources and sinks") + + p2E = plot(μ * t, [dEdt_computed, dEdt_numerical], + label = ["computed Wᵋ-Dᵋ" "numerical dE/dt"], + linestyle = [:solid :dashdotdot], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "dE/dt") + + p3E = plot(μ * t, residual_E, + label = "residual dE/dt = computed - numerical", + linewidth = 2, + alpha = 0.7, + xlabel = "μt") + + t = Z.t[2:E.i] + + p1Z = plot(μ * t, [Wᶻ[2:Z.i] εᶻ.+0*t -Dᶻ[1:Z.i-1] -Rᶻ[1:Z.i-1]], + label = ["enstrophy work, Wᶻ" "mean enstrophy work, " "enstrophy dissipation, Dᶻ" "enstrophy drag, Rᶻ = - 2μZ"], + linestyle = [:solid :dash :solid :solid], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "enstrophy sources and sinks") + + + p2Z = plot(μ * t, [dZdt_computed, dZdt_numerical], + label = ["computed Wᶻ-Dᶻ" "numerical dZ/dt"], + linestyle = [:solid :dashdotdot], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "dZ/dt") + + p3Z = plot(μ * t, residual_Z, + label = "residual dZ/dt = computed - numerical", + linewidth = 2, + alpha = 0.7, + xlabel = "μt") + + layout = @layout Plots.grid(3, 2) + + pbudgets = plot(p1E, p1Z, p2E, p2Z, p3E, p3Z, layout=layout, size = (900, 1200)) + + return pζ, pbudgets +end +nothing # hide + + + + +# ## Time-stepping the `Problem` forward + +# Finally, we time-step the `Problem` forward in time. + +startwalltime = time() +for i = 1:ns + stepforward!(prob, diags, round(Int, nt/ns)) + + TwoDNavierStokes.updatevars!(prob) + + cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) + + log = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min", clock.step, clock.t, + cfl, (time()-startwalltime)/60) + + println(log) +end + + +# ## Plot +# And now let's see what we got. We plot the output. First the final snapshot +# of the vorticity field. + +pζ, pbudgets = computetendencies_and_makeplot(prob, diags) + +pζ + +# And finaly the energy and enstrophy budgets. + +pbudgets From 0a51dffa8cc19d64d40e999c96ae6313b5f5b43f Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 09:53:10 +1000 Subject: [PATCH 3/9] adds a 2DNS forced example with an animation --- docs/make.jl | 2 + .../twodnavierstokes_stochasticforcing.jl | 173 ++++++++++++++++++ ...dnavierstokes_stochasticforcing_budgets.jl | 2 +- 3 files changed, 176 insertions(+), 1 deletion(-) create mode 100644 examples/twodnavierstokes_stochasticforcing.jl diff --git a/docs/make.jl b/docs/make.jl index a96d1f4c..f9073745 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -20,6 +20,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", @@ -72,6 +73,7 @@ sitename = "GeophysicalFlows.jl", "Examples" => [ "TwoDNavierStokes" => Any[ "generated/twodnavierstokes_decaying.md", + "generated/twodnavierstokes_stochasticforcing.md", "generated/twodnavierstokes_stochasticforcing_budgets.md", ], "BarotropicQG" => Any[ diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl new file mode 100644 index 00000000..81d3cb75 --- /dev/null +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -0,0 +1,173 @@ +# # 2D forced-dissipative turbulence +# +#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +# +# A simulation of forced-dissipative two-dimensional turbulence. We solve the +# two-dimensional vorticity equation with stochastic excitation and dissipation in +# the form of linear drag and hyperviscosity. + +using FourierFlows, Printf, Plots + +using FourierFlows: parsevalsum +using Random: seed! +using FFTW: irfft + +import GeophysicalFlows.TwoDNavierStokes +import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy + + +# ## Choosing a device: CPU or GPU + +dev = CPU() # Device (CPU/GPU) +nothing # hide + + +# ## Numerical, domain, and simulation parameters +# +# First, we pick some numerical and physical parameters for our model. + + n, L = 256, 2π # grid resolution and domain length + ν, nν = 2e-7, 2 # hyperviscosity coefficient and order + μ, nμ = 1e-1, 0 # linear drag coefficient + dt = 0.005 # timestep and final time +nsteps = 4000 # total number of steps + nsubs = 20 # number of steps between each plot +nothing # hide + + +# ## Forcing + +# We force the vorticity equation with stochastic excitation that is delta-correlated +# in time and while spatially homogeneously and isotropically correlated. The forcing +# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and +# width $\delta k_f$, and it injects energy per unit area and per unit time equal +# to $\varepsilon$. + +forcing_wavenumber = 14.0 # the central forcing wavenumber for a spectrum that is a ring in wavenumber space +forcing_bandwidth = 1.5 # the width of the forcing spectrum +ε = 0.1 # energy input rate by the forcing + +grid = TwoDGrid(dev, n, L) +x, y = grid.x, grid.y + +K = sqrt.(grid.Krsq) +forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2)) +forcing_spectrum[K .< (2*2π/L)] .= 0 # make sure that focing has no power for low wavenumbers +forcing_spectrum[K .> (20*2π/L)] .= 0 # make sure that focing has no power for high wavenumbers +ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) +@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy ε + +seed!(1234) +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)) + ξ[1, 1] = 0 + @. Fh = ξ * sqrt(forcing_spectrum) + nothing +end +nothing # hide + + +# ## Problem setup +# We initialize a `Problem` by providing a set of keyword arguments. The +# `stepper` keyword defines the time-stepper to be used. +prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4", + calcF=calcF!, stochastic=true) +nothing # hide + +# Define some shortcuts for convenience. +sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid +nothing # hide + + +# First let's see how a forcing realization looks like. +calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) + +heatmap(x, y, irfft(vars.Fh, grid.nx), + aspectratio = 1, + c = :balance, + clim = (-200, 200), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "x", + ylabel = "y", + title = "a forcing realization", + framestyle = :box) + + +# ## Setting initial conditions + +# Our initial condition is simply fluid at rest. +TwoDNavierStokes.set_zeta!(prob, zeros(grid.nx, grid.ny)) + + +# ## Diagnostics + +# Create Diagnostics; the diagnostics are aimed to probe the energy budget. +E = Diagnostic(energy, prob, nsteps=nsteps) # energy +Z = Diagnostic(enstrophy, prob, nsteps=nsteps) # enstrophy +diags = [E, Z] # a list of Diagnostics passed to `stepforward!` will be updated every timestep. +nothing # hide + + +# ## Visualizing the simulation + +# We initialize a plot with the vorticity field and the time-series of +# energy and enstrophy diagnostics. To plot energy and enstrophy on the same +# axes we scale enstrophy with $k_f^2$. + +p1 = heatmap(x, y, vars.zeta, + aspectratio = 1, + c = :balance, + clim = (-40, 40), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "x", + ylabel = "y", + title = "vorticity, t="*@sprintf("%.2f", clock.t), + framestyle = :box) + +p2 = plot(2, # this means "a plot with two series" + label = ["energy E(t)" "enstrophy Z(t) / k_f²"], + legend = :right, + linewidth = 2, + alpha = 0.7, + xlabel = "μ t", + xlims = (0, 1.1 * μ * nsteps * dt), + ylims = (0, 0.6)) + +l = @layout Plots.grid(1, 2) +p = plot(p1, p2, layout = l, size = (900, 420)) + + +# ## Time-stepping the `Problem` forward + +# Finally, we time-step the `Problem` forward in time. + +startwalltime = time() +anim = @animate for j = 0:Int(nsteps/nsubs) + + log = @sprintf("step: %04d, t: %d, E: %.4f, Z: %.4f, walltime: %.2f min", + clock.step, clock.t, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60) + + if j%(1000/nsubs)==0; println(log) end + + p[1][1][:z] = vars.zeta + p[1][:title] = "vorticity, μ t = "*@sprintf("%.2f", clock.t) + push!(p[2][1], μ * E.t[E.i], E.data[E.i]) + push!(p[2][2], μ * Z.t[Z.i], Z.data[Z.i] / forcing_wavenumber^2) + + stepforward!(prob, diags, nsubs) + TwoDNavierStokes.updatevars!(prob) + +end + +mp4(anim, "twodturb_forced.mp4", fps=18) +gif(anim, "twodturb_forced.gif", fps=18) \ No newline at end of file diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index ce832f63..936b298d 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -1,4 +1,4 @@ -# # 2D forced-dissipative turbulence +# # 2D forced-dissipative turbulence budgets # #md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). #md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). From e3efe3187ef14be18d4f088d5eed8ec3685dad0b Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 10:18:05 +1000 Subject: [PATCH 4/9] minor clean up of the forcing spectrum construction --- examples/barotropicqg_betaforced.jl | 17 +++++++++-------- examples/barotropicqgql_betaforced.jl | 17 +++++++++-------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/examples/barotropicqg_betaforced.jl b/examples/barotropicqg_betaforced.jl index 1a711939..981a2d2d 100644 --- a/examples/barotropicqg_betaforced.jl +++ b/examples/barotropicqg_betaforced.jl @@ -55,14 +55,15 @@ forcing_bandwidth = 1.5 # the width of the forcing spectrum gr = TwoDGrid(nx, Lx) -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π/Lx ] .= 0 # make sure forcing does not have power at k=0 +ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) +@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy at rate ε seed!(1234) # reset of the random number generator for reproducibility nothing # hide @@ -72,7 +73,7 @@ function calcFq!(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 + return nothing end nothing # hide diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index e306042e..0f002297 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -56,14 +56,15 @@ forcing_bandwidth = 1.5 # the width of the forcing spectrum gr = TwoDGrid(nx, Lx) -k = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl ] # a 2D grid 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 +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(-(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π/Lx ] .= 0 # make sure forcing does not have power at k=0 +ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) +@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy at rate ε seed!(1234) # reset of the random number generator for reproducibility nothing # hide From 9c7449b1707c41e418e687ab74242c7ac464a0b7 Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 10:30:16 +1000 Subject: [PATCH 5/9] some clarifications and comments --- .../twodnavierstokes_stochasticforcing.jl | 25 ++++++------ ...dnavierstokes_stochasticforcing_budgets.jl | 39 ++++++++++--------- 2 files changed, 35 insertions(+), 29 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 81d3cb75..cbdd1ca4 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -28,11 +28,11 @@ nothing # hide # First, we pick some numerical and physical parameters for our model. n, L = 256, 2π # grid resolution and domain length - ν, nν = 2e-7, 2 # hyperviscosity coefficient and order + ν, nν = 2e-7, 2 # hyperviscosity coefficient and hyperviscosity order μ, nμ = 1e-1, 0 # linear drag coefficient - dt = 0.005 # timestep and final time -nsteps = 4000 # total number of steps - nsubs = 20 # number of steps between each plot + dt = 0.005 # timestep +nsteps = 4000 # total number of steps + nsubs = 20 # number of steps between each plot nothing # hide @@ -49,14 +49,13 @@ forcing_bandwidth = 1.5 # the width of the forcing spectrum ε = 0.1 # energy input rate by the forcing grid = TwoDGrid(dev, n, L) -x, y = grid.x, grid.y -K = sqrt.(grid.Krsq) +K = @. sqrt(grid.Krsq) forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2)) -forcing_spectrum[K .< (2*2π/L)] .= 0 # make sure that focing has no power for low wavenumbers -forcing_spectrum[K .> (20*2π/L)] .= 0 # make sure that focing has no power for high wavenumbers +forcing_spectrum[K .< ( 2 * 2π/L)] .= 0 # no power at low wavenumbers +forcing_spectrum[K .> (20 * 2π/L)] .= 0 # no power at high wavenumbers ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) -@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy ε +@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy at rate ε seed!(1234) nothing # hide @@ -80,10 +79,14 @@ nothing # hide # Define some shortcuts for convenience. sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + +x, y = grid.x, grid.y nothing # hide -# First let's see how a forcing realization looks like. +# First let's see how a forcing realization looks like. Function `calcF!()` computes +# the forcing in Fourier space and saves it into variable `vars.Fh`, so we first need to +# go back to physical space. calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) heatmap(x, y, irfft(vars.Fh, grid.nx), @@ -102,7 +105,7 @@ heatmap(x, y, irfft(vars.Fh, grid.nx), # ## Setting initial conditions -# Our initial condition is simply fluid at rest. +# Our initial condition is a fluid at rest. TwoDNavierStokes.set_zeta!(prob, zeros(grid.nx, grid.ny)) diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 936b298d..3d6370f6 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -30,12 +30,12 @@ nothing # hide # # First, we pick some numerical and physical parameters for our model. - n, L = 256, 2π # grid resolution and domain length - ν, nν = 2e-7, 2 # hyperviscosity coefficient and order - μ, nμ = 1e-1, 0 # linear drag coefficient -dt, tf = 0.005, 0.2 / μ # timestep and final time - nt = round(Int, tf / dt) # total timesteps - ns = 4 # how many intermediate times we want to plot + n, L = 256, 2π # grid resolution and domain length + ν, nν = 2e-7, 2 # hyperviscosity coefficient and hyperviscosity order + μ, nμ = 1e-1, 0 # linear drag coefficient +dt, tf = 0.005, 0.2 / μ # timestep and final time + nt = round(Int, tf / dt) # total timesteps + ns = 4 # how many intermediate times we want to plot nothing # hide @@ -52,14 +52,13 @@ forcing_bandwidth = 1.5 # the width of the forcing spectrum ε = 0.1 # energy input rate by the forcing grid = TwoDGrid(dev, n, L) -x, y = grid.x, grid.y -K = sqrt.(grid.Krsq) +K = @. sqrt(grid.Krsq) forcing_spectrum = @. exp(-(K - forcing_wavenumber)^2 / (2 * forcing_bandwidth^2)) -forcing_spectrum[K .< (2*2π/L)] .= 0 # make sure that focing has no power for low wavenumbers -forcing_spectrum[K .> (20*2π/L)] .= 0 # make sure that focing has no power for high wavenumbers +forcing_spectrum[K .< ( 2 * 2π/L)] .= 0 # no power at low wavenumbers +forcing_spectrum[K .> (20 * 2π/L)] .= 0 # no power at high wavenumbers ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) -@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy ε +@. forcing_spectrum *= ε/ε0 # normalize forcing to inject energy at rate ε seed!(1234) nothing # hide @@ -83,10 +82,14 @@ nothing # hide # Define some shortcuts for convenience. sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + +x, y = grid.x, grid.y nothing # hide -# First let's see how a forcing realization looks like. +# First let's see how a forcing realization looks like. Function `calcF!()` computes +# the forcing in Fourier space and saves it into variable `vars.Fh`, so we first need to +# go back to physical space. calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) heatmap(x, y, irfft(vars.Fh, grid.nx), @@ -105,13 +108,13 @@ heatmap(x, y, irfft(vars.Fh, grid.nx), # ## Setting initial conditions -# Our initial condition is simply fluid at rest. +# Our initial condition is a fluid at rest. TwoDNavierStokes.set_zeta!(prob, zeros(grid.nx, grid.ny)) # ## Diagnostics -# Create Diagnostics; the diagnostics are aimed to probe the energy budget. +# Create Diagnostics; the diagnostics are aimed to probe the energy and enstrophy budgets. E = Diagnostic(energy, prob, nsteps=nt) # energy Rᵋ = Diagnostic(energy_drag, prob, nsteps=nt) # energy dissipation by drag Dᵋ = Diagnostic(energy_dissipation, prob, nsteps=nt) # energy dissipation by hyperviscosity @@ -127,9 +130,9 @@ nothing # hide # ## Visualizing the simulation # We define a function that plots the vorticity field and the evolution of -# the diagnostics: energy, enstrophy and all terms involved in the energy and -# enstrophy budgets. Last, we confirm whether the energy and enstrophy budgets -# are accurately computed, e.g., $\mathrm{d}E/\mathrm{d}t = W^\varepsilon - +# the diagnostics: energy, enstrophy, and all terms involved in the energy and +# enstrophy budgets. Last, we also check (by plotting) whether the energy and enstrophy +# budgets are accurately computed, e.g., $\mathrm{d}E/\mathrm{d}t = W^\varepsilon - # R^\varepsilon - D^\varepsilon$. function computetendencies_and_makeplot(prob, diags) @@ -248,7 +251,7 @@ end # ## Plot -# And now let's see what we got. We plot the output. First the final snapshot +# And now let's see what we got. First we plot the final snapshot # of the vorticity field. pζ, pbudgets = computetendencies_and_makeplot(prob, diags) From de0c2e961dbbee084fadfd6ab59003392182a0ef Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 10:49:12 +1000 Subject: [PATCH 6/9] small fixes --- examples/barotropicqg_betaforced.jl | 18 +++++++++--------- examples/barotropicqgql_betaforced.jl | 24 ++++++++++++------------ 2 files changed, 21 insertions(+), 21 deletions(-) diff --git a/examples/barotropicqg_betaforced.jl b/examples/barotropicqg_betaforced.jl index 981a2d2d..0efc9806 100644 --- a/examples/barotropicqg_betaforced.jl +++ b/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 @@ -53,16 +53,16 @@ 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 = @. 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(-(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π/Lx ] .= 0 # make sure forcing does not have power at k=0 -ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) +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 @@ -70,9 +70,9 @@ 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 + Fh[abs.(grid.Krsq) .== 0] .= 0 return nothing end nothing # hide @@ -83,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 diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 0f002297..809e8bf8 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/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,7 +36,7 @@ nothing # hide # ## Physical parameters -Lx = 2π # domain size + L = 2π # domain size β = 10.0 # planetary PV gradient μ = 0.01 # bottom drag nothing # hide @@ -54,16 +54,16 @@ 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 = @. 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(-(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π/Lx ] .= 0 # make sure forcing does not have power at k=0 -ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly) +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 @@ -71,10 +71,10 @@ 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 @@ -84,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 From d21a6a6d772a05bd27de01b7f28c5bdcd7f7ee8b Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 11:02:39 +1000 Subject: [PATCH 7/9] fix nbviewer/binder urls + stop creating gif --- examples/twodnavierstokes_stochasticforcing.jl | 1 - examples/twodnavierstokes_stochasticforcing_budgets.jl | 4 ++-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index cbdd1ca4..cc832746 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -173,4 +173,3 @@ anim = @animate for j = 0:Int(nsteps/nsubs) end mp4(anim, "twodturb_forced.mp4", fps=18) -gif(anim, "twodturb_forced.gif", fps=18) \ No newline at end of file diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 3d6370f6..a4a630ff 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -1,7 +1,7 @@ # # 2D forced-dissipative turbulence budgets # -#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). -#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing_budgets.ipynb). +#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing_budgets.ipynb). # # A simulation of forced-dissipative two-dimensional turbulence. We solve the # two-dimensional vorticity equation with stochastic excitation and dissipation in From a8be81845cc671d4ef414ce64d34539cce000b3a Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 11:32:07 +1000 Subject: [PATCH 8/9] still there were some typos... --- examples/barotropicqg_betadecay.jl | 25 +++++++++++++------------ examples/barotropicqg_betaforced.jl | 12 ++++++------ examples/barotropicqgql_betaforced.jl | 6 +++--- examples/multilayerqg_2layer.jl | 6 +++--- 4 files changed, 25 insertions(+), 24 deletions(-) diff --git a/examples/barotropicqg_betadecay.jl b/examples/barotropicqg_betadecay.jl index ed7fd2b0..851dc244 100644 --- a/examples/barotropicqg_betadecay.jl +++ b/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,9 +32,9 @@ 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 @@ -42,7 +42,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) 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) diff --git a/examples/barotropicqg_betaforced.jl b/examples/barotropicqg_betaforced.jl index 0efc9806..61069419 100644 --- a/examples/barotropicqg_betaforced.jl +++ b/examples/barotropicqg_betaforced.jl @@ -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 @@ -55,13 +55,13 @@ forcing_bandwidth = 1.5 # the width of the forcing spectrum gr = TwoDGrid(n, L) -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 +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(-(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 +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 ε diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 809e8bf8..265bc92c 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -36,9 +36,9 @@ nothing # hide # ## Physical parameters - L = 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 diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 8f9238d1..afccef5c 100644 --- a/examples/multilayerqg_2layer.jl +++ b/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 From 07b3b61092c60c3c66c61f175b513fcbd8c188d4 Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Sat, 12 Sep 2020 13:06:54 +1000 Subject: [PATCH 9/9] fix units in plot title --- examples/twodnavierstokes_stochasticforcing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index cc832746..1995b431 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -144,7 +144,7 @@ p2 = plot(2, # this means "a plot with two series" alpha = 0.7, xlabel = "μ t", xlims = (0, 1.1 * μ * nsteps * dt), - ylims = (0, 0.6)) + ylims = (0, 0.55)) l = @layout Plots.grid(1, 2) p = plot(p1, p2, layout = l, size = (900, 420)) @@ -163,7 +163,7 @@ anim = @animate for j = 0:Int(nsteps/nsubs) if j%(1000/nsubs)==0; println(log) end p[1][1][:z] = vars.zeta - p[1][:title] = "vorticity, μ t = "*@sprintf("%.2f", clock.t) + p[1][:title] = "vorticity, μ t = "*@sprintf("%.2f", μ * clock.t) push!(p[2][1], μ * E.t[E.i], E.data[E.i]) push!(p[2][2], μ * Z.t[Z.i], Z.data[Z.i] / forcing_wavenumber^2)