diff --git a/Project.toml b/Project.toml index d8e56b02..f0b6ff2e 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" CUDA = "^1, ^2.4.2, 3.0.0 - 3.6.4, ^3.7.1" DocStringExtensions = "^0.8, 0.9" FFTW = "^1" -FourierFlows = "^0.9.2" +FourierFlows = "^0.9.4" JLD2 = "^0.1, ^0.2, ^0.3, ^0.4" Reexport = "^0.2, ^1" SpecialFunctions = "^0.10, ^1, 2" diff --git a/README.md b/README.md index 7b211aac..e8d4ced5 100644 --- a/README.md +++ b/README.md @@ -62,14 +62,17 @@ Some animations created with GeophysicalFlows.jl are [online @ youtube]. ## Modules -* `TwoDNavierStokes`: the two-dimensional vorticity equation. -* `SingleLayerQG`: the barotropic or equivalent-barotropic quasi-geostrophic equation, which - generalizes `TwoDNavierStokes` to cases with topography, Coriolis parameters of the form - `f = f₀ + βy`, and finite Rossby radius of deformation. -* `MultiLayerQG`: a multi-layer quasi-geostrophic model over topography and with the ability - to impose a zonal flow `U_n(y)` in each layer. -* `SurfaceQG`: a surface quasi-geostrophic model. -* `BarotropicQGQL`: the quasi-linear barotropic quasi-geostrophic equation. +* [`TwoDNavierStokes`](https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/modules/twodnavierstokes/): the + two-dimensional vorticity equation. +* [`SingleLayerQG`](https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/modules/singlelayerqg/): the barotropic + or equivalent-barotropic quasi-geostrophic equation, which generalizes `TwoDNavierStokes` to cases with topography, Coriolis + parameters of the form `f = f₀ + βy`, and finite Rossby radius of deformation. +* [`MultiLayerQG`](https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/modules/multilayerqg/): a multi-layer + quasi-geostrophic model over topography allowing to impose a zonal flow `U_n(y)` in each layer. +* [`SurfaceQG`](https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/modules/surfaceqg/): a surface + quasi-geostrophic model. +* [`BarotropicQGQL`](https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/modules/barotropicqgql/): the + quasi-linear barotropic quasi-geostrophic equation. ## Scalability @@ -143,7 +146,6 @@ than happy to help along the way. For more information, check out our [contributors' guide](https://github.com/FourierFlows/GeophysicalFlows.jl/blob/main/CONTRIBUTING.md). - [FourierFlows.jl]: https://github.com/FourierFlows/FourierFlows.jl [documentation]: https://fourierflows.github.io/GeophysicalFlowsDocumentation/dev/ [online @ youtube]: https://www.youtube.com/channel/UCO_0ugkNUwCsFUMtepwYTqw diff --git a/docs/Project.toml b/docs/Project.toml index 0511b5cd..e0a67d1e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,11 +1,11 @@ [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" [compat] CUDA = "^1, ^2.4.2, 3.0.0 - 3.6.4, ^3.7.1" Literate = "≥2.9.0" -Plots = "1.10.1 - 1.21.3" diff --git a/docs/make.jl b/docs/make.jl index b4148021..dcc5782a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,14 +1,9 @@ using Documenter, Literate, - Plots, # to not capture precompilation output + CairoMakie, # to not capture precompilation output GeophysicalFlows -# Gotta set this environment variable when using the GR run-time on CI machines. -# This happens as examples will use Plots.jl to make plots and movies. -# See: https://github.com/jheinen/GR.jl/issues/278 -ENV["GKSwstype"] = "100" - ##### ##### Generate literated examples ##### diff --git a/docs/src/modules/barotropicqgql.md b/docs/src/modules/barotropicqgql.md index 7db16a6c..ca08d87e 100644 --- a/docs/src/modules/barotropicqgql.md +++ b/docs/src/modules/barotropicqgql.md @@ -115,4 +115,4 @@ Other diagnostic include: [`dissipation`](@ref GeophysicalFlows.BarotropicQGQL.d ## Examples -- [`examples/barotropicqgql_betaforced.jl`](../../literated/barotropicqgql_betaforced/): A script that simulates forced-dissipative quasi-linear quasi-geostrophic flow on a beta plane demonstrating zonation. The forcing is temporally delta-correlated and its spatial structure is isotropic with power in a narrow annulus of total radius ``k_f`` in wavenumber space. This example demonstrates that the anisotropic inverse energy cascade is not required for zonation. +- [`examples/barotropicqgql_betaforced.jl`](@ref barotropicqgql_betaforced_example): Simulate forced-dissipative quasi-linear quasi-geostrophic flow on a beta plane demonstrating zonation. The forcing is temporally delta-correlated and its spatial structure is isotropic with power in a narrow annulus of total radius ``k_f`` in wavenumber space. This example demonstrates that the anisotropic inverse energy cascade is not required for zonation. diff --git a/docs/src/modules/multilayerqg.md b/docs/src/modules/multilayerqg.md index 650e0d3c..313736c4 100644 --- a/docs/src/modules/multilayerqg.md +++ b/docs/src/modules/multilayerqg.md @@ -146,4 +146,4 @@ GeophysicalFlows.MultiLayerQG.fluxes ## Examples - - [`examples/multilayerqg_2layer.jl`](../../literated/multilayerqg_2layer/): A script that simulates the growth and equilibration of baroclinic eddy turbulence in the Phillips 2-layer model. + - [`examples/multilayerqg_2layer.jl`](@ref multilayerqg_2layer_example): Simulate the growth and equilibration of baroclinic eddy turbulence in the Phillips 2-layer model. diff --git a/docs/src/modules/singlelayerqg.md b/docs/src/modules/singlelayerqg.md index cd3a0533..6590d839 100644 --- a/docs/src/modules/singlelayerqg.md +++ b/docs/src/modules/singlelayerqg.md @@ -110,10 +110,10 @@ Other diagnostic include: [`energy_dissipation`](@ref GeophysicalFlows.SingleLay ## Examples -- [`examples/singlelayerqg_betadecay.jl`](../../literated/singlelayerqg_betadecay/): A script that simulates decaying quasi-geostrophic flow on a beta plane demonstrating zonation. +- [`examples/singlelayerqg_betadecay.jl`](@ref singlelayerqg_betadecay_example): Simulate decaying quasi-geostrophic flow on a beta plane demonstrating zonation. -- [`examples/singlelayerqg_betaforced.jl`](../../literated/singlelayerqg_betaforced/): A script that simulates forced-dissipative quasi-geostrophic flow on a beta plane demonstrating zonation. The forcing is temporally delta-correlated with isotropic spatial structure with power in a narrow annulus in wavenumber space with total wavenumber ``k_f``. +- [`examples/singlelayerqg_betaforced.jl`](@ref singlelayerqg_betaforced_example): Simulate forced-dissipative quasi-geostrophic flow on a beta plane demonstrating zonation. The forcing is temporally delta-correlated with isotropic spatial structure with power in a narrow annulus in wavenumber space with total wavenumber ``k_f``. -- [`examples/singlelayerqg_decay_topography.jl`](../../literated/singlelayerqg_decay_topography/): A script that simulates two dimensional turbulence (barotropic quasi-geostrophic flow with ``\beta=0``) above topography. +- [`examples/singlelayerqg_decay_topography.jl`](@ref singlelayerqg_decay_topography_example): Simulate two dimensional turbulence (barotropic quasi-geostrophic flow with ``\beta=0``) above topography. -- [`examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl`](../../literated singlelayerqg_decaying_barotropic_equivalentbarotropic/): A script that simulates two dimensional turbulence (``\beta=0``) with both infinite and finite Rossby radius of deformation and compares the evolution of the two. \ No newline at end of file +- [`examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl`](@ref singlelayerqg_decaying_barotropic_equivalentbarotropic_example): Simulate two dimensional turbulence (``\beta=0``) with both infinite and finite Rossby radius of deformation and compares the evolution of the two. \ No newline at end of file diff --git a/docs/src/modules/surfaceqg.md b/docs/src/modules/surfaceqg.md index 552c95ba..f40b3f67 100644 --- a/docs/src/modules/surfaceqg.md +++ b/docs/src/modules/surfaceqg.md @@ -96,6 +96,6 @@ Other diagnostic include: [`buoyancy_dissipation`](@ref GeophysicalFlows.Surface ## Examples -- [`examples/surfaceqg_decaying.jl`](../../literated/surfaceqg_decaying/): A script that simulates decaying surface quasi-geostrophic flow with a prescribed initial buoyancy field, producing an animation of the evolution of the surface buoyancy. +- [`examples/surfaceqg_decaying.jl`](@ref surfaceqg_decaying_example): Simulate decaying surface quasi-geostrophic flow with a prescribed initial buoyancy field. > Capet, X. et al., (2008). Surface kinetic energy transfer in surface quasi-geostrophic flows. *J. Fluid Mech.*, **604**, 165-174. diff --git a/docs/src/modules/twodnavierstokes.md b/docs/src/modules/twodnavierstokes.md index a678a4e3..23a65ffc 100644 --- a/docs/src/modules/twodnavierstokes.md +++ b/docs/src/modules/twodnavierstokes.md @@ -90,10 +90,10 @@ Other diagnostic include: [`energy_dissipation`](@ref GeophysicalFlows.TwoDNavie ## Examples -- [`examples/twodnavierstokes_decaying.jl`](../../literated/twodnavierstokes_decaying/): A script that simulates decaying two-dimensional turbulence reproducing the results by +- [`examples/twodnavierstokes_decaying.jl`](@ref twodnavierstokes_decaying_example): Simulates decaying two-dimensional turbulence reproducing the results by: > McWilliams, J. C. (1984). The emergence of isolated coherent vortices in turbulent flow. *J. Fluid Mech.*, **146**, 21-43. -- [`examples/twodnavierstokes_stochasticforcing.jl`](../../literated/twodnavierstokes_stochasticforcing/): A script that simulates forced-dissipative two-dimensional turbulence with isotropic temporally delta-correlated stochastic forcing. +- [`examples/twodnavierstokes_stochasticforcing.jl`](@ref twodnavierstokes_stochasticforcing_example): Simulate forced-dissipative two-dimensional turbulence with isotropic temporally delta-correlated stochastic forcing. -- [`examples/twodnavierstokes_stochasticforcing_budgets.jl`](../../literated/twodnavierstokes_stochasticforcing_budgets/): A script that simulates forced-dissipative two-dimensional turbulence demonstrating how we can compute the energy and enstrophy budgets. +- [`examples/twodnavierstokes_stochasticforcing_budgets.jl`](@ref twodnavierstokes_stochasticforcing_budgets_example): Simulate forced-dissipative two-dimensional turbulence demonstrating how we can compute the energy and enstrophy budgets. diff --git a/examples/Project.toml b/examples/Project.toml index 1669ce76..28c4c8e3 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -1,11 +1,8 @@ [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" - -[compat] -Plots = "≥ 1.10.1" diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 8259fbba..23765e51 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -1,4 +1,4 @@ -# # Quasi-Linear forced-dissipative barotropic QG beta-plane turbulence +# # [Quasi-Linear forced-dissipative barotropic QG beta-plane turbulence](@id barotropicqgql_betaforced_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/barotropicqgql_betaforced.ipynb). # @@ -12,15 +12,18 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, CUDA, Random, Statistics, Printf, Plots" +# pkg"add GeophysicalFlows, CUDA, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. + +using GeophysicalFlows, CUDA, Random, Printf, CairoMakie -using GeophysicalFlows, CUDA, Random, Printf, Plots using Statistics: mean + parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU @@ -95,33 +98,36 @@ nothing # hide # at every time-step that removes enstrophy at high wavenumbers and, thereby, # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. # Thus, we choose not to do any dealiasing by providing `aliased_fraction=0`. -prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper, +prob = BarotropicQGQL.Problem(dev; nx=n, Lx=L, β, μ, dt, stepper, calcF=calcF!, stochastic=true, aliased_fraction=0) nothing # hide # and define some shortcuts. sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -x, y = grid.x, grid.y +x, y = grid.x, grid.y +Lx, Ly = grid.Lx, grid.Ly nothing # hide # First let's see how a forcing realization looks like. Note that when plotting, we decorate # the variable to be plotted with `Array()` to make sure it is brought back on the CPU when # `vars` live on the GPU. + calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) -heatmap(x, y, Array(irfft(vars.Fh, grid.nx)'), - aspectratio = 1, - c = :balance, - clim = (-8, 8), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, +fig = Figure() + +ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", - title = "a forcing realization", - framestyle = :box) + aspect = 1, + title = "a forcing realization", + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); + colormap = :balance, colorrange = (-8, 8)) + +fig # ## Setting initial conditions @@ -133,17 +139,14 @@ nothing # hide # ## Diagnostics # Create Diagnostics -- `energy` and `enstrophy` are functions imported at the top. -E = Diagnostic(BarotropicQGQL.energy, prob; nsteps=nsteps) -Z = Diagnostic(BarotropicQGQL.enstrophy, prob; nsteps=nsteps) +E = Diagnostic(BarotropicQGQL.energy, prob; nsteps) +Z = Diagnostic(BarotropicQGQL.enstrophy, prob; nsteps) nothing # hide # We can also define our custom diagnostics via functions. -function zetaMean(prob) - sol = prob.sol - sol[1, :] -end +zetaMean(prob) = prob.sol[1, :] -zMean = Diagnostic(zetaMean, prob; nsteps=nsteps, freq=10) # the zonal-mean vorticity +zMean = Diagnostic(zetaMean, prob; nsteps, freq=10) # the zonal-mean vorticity nothing # hide # We combile all diags in a list. @@ -166,8 +169,17 @@ if !isdir(plotpath); mkdir(plotpath); end nothing # hide # and then create Output. -get_sol(prob) = sol # extracts the Fourier-transformed solution -get_u(prob) = irfft(im * g.l .* g.invKrsq .* sol, g.nx) +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution + +function get_u(prob) + grid, vars = prob.grid, prob.vars + + @. vars.uh = im * grid.l * grid.invKrsq * sol + ldiv!(vars.u, grid.rfftplan, deepcopy(vars.uh)) + + return vars.u +end + out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) @@ -177,99 +189,82 @@ out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) # corresponding zonal-mean vorticity and zonal-mean zonal velocity and timeseries # of energy and enstrophy. -function plot_output(prob) - ζ̄, ζ′= prob.vars.Zeta, prob.vars.zeta - ζ = @. ζ̄ + ζ′ - ψ̄, ψ′= prob.vars.Psi, prob.vars.psi - ψ = @. ψ̄ + ψ′ - ζ̄ₘ = mean(ζ̄, dims=1)' - ūₘ = mean(prob.vars.U, dims=1)' - - pζ = heatmap(x, y, Array(ζ'), - aspectratio = 1, - legend = false, - c = :balance, - clim = (-8, 8), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "vorticity ζ=∂v/∂x-∂u/∂y", - framestyle = :box) - - pψ = contourf(x, y, Array(ψ'), - levels = -0.32:0.04:0.32, - aspectratio = 1, - linewidth = 1, - legend = false, - clim = (-0.22, 0.22), - c = :viridis, - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "streamfunction ψ", - framestyle = :box) - - pζm = plot(Array(ζ̄ₘ), y, - legend = false, - linewidth = 2, - alpha = 0.7, - yticks = -3:3, - xlims = (-3, 3), - xlabel = "zonal mean ζ", - ylabel = "y") - plot!(pζm, 0*y, y, linestyle=:dash, linecolor=:black) - - pum = plot(Array(ūₘ), y, - legend = false, - linewidth = 2, - alpha = 0.7, - yticks = -3:3, - xlims = (-0.5, 0.5), - xlabel = "zonal mean u", - ylabel = "y") - plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black) - - pE = plot(1, - label = "energy", - linewidth = 2, - alpha = 0.7, - xlims = (-0.1, 4.1), - ylims = (0, 0.05), - xlabel = "μt") - - pZ = plot(1, - label = "enstrophy", - linecolor = :red, - legend = :bottomright, - linewidth = 2, - alpha = 0.7, - xlims = (-0.1, 4.1), - ylims = (0, 5), - xlabel = "μt") - - l = @layout Plots.grid(2, 3) - p = plot(pζ, pζm, pE, pψ, pum, pZ, layout=l, size = (1000, 600)) - - return p -end +title_ζ = Observable(@sprintf("vorticity, μt = %.2f", μ * clock.t)) +title_ψ = "streamfunction ψ" + +fig = Figure(resolution=(1000, 600)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axζ = Axis(fig[1, 1]; title = title_ζ, axis_kwargs...) + +axψ = Axis(fig[2, 1]; title = title_ψ, axis_kwargs...) + +axζ̄ = Axis(fig[1, 2], + xlabel = "zonal mean ζ", + ylabel = "y", + aspect = 1, + limits = ((-3, 3), (-Ly/2, Ly/2))) + +axū = Axis(fig[2, 2], + xlabel = "zonal mean u", + ylabel = "y", + aspect = 1, + limits = ((-0.5, 0.5), (-Ly/2, Ly/2))) + +axE = Axis(fig[1, 3], + xlabel = "μ t", + ylabel = "energy", + aspect = 1, + limits = ((-0.1, 4.1), (0, 0.05))) + +axZ = Axis(fig[2, 3], + xlabel = "μ t", + ylabel = "enstrophy", + aspect = 1, + limits = ((-0.1, 4.1), (0, 5))) + +ζ̄, ζ′= prob.vars.Zeta, prob.vars.zeta +ζ = Observable(Array(@. ζ̄ + ζ′)) +ψ̄, ψ′= prob.vars.Psi, prob.vars.psi +ψ = Observable(Array(@. ψ̄ + ψ′)) +ζ̄ₘ = Observable(Array(vec(mean(ζ̄, dims=1)))) +ūₘ = Observable(Array(vec(mean(prob.vars.U, dims=1)))) + +μt = Observable(μ * E.t[1:1]) +energy = Observable(E.data[1:1]) +enstrophy = Observable(Z.data[1:1]) + +heatmap!(axζ, x, y, ζ; + colormap = :balance, colorrange = (-8, 8)) + +heatmap!(axψ, x, y, ψ; + colormap = :viridis, colorrange = (-0.22, 0.22)) + +lines!(axζ̄, ζ̄ₘ, y; linewidth = 3) +lines!(axζ̄, 0y, y; linewidth = 1, linestyle=:dash) + +lines!(axū, ūₘ, y; linewidth = 3) +lines!(axū, 0y, y; linewidth = 1, linestyle=:dash) + +lines!(axE, μt, energy; linewidth = 3) +lines!(axZ, μt, enstrophy; linewidth = 3, color = :red) + nothing # hide # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. - -p = plot_output(prob) +# We step the `Problem` forward in time. startwalltime = time() -anim = @animate for j = 0:round(Int, nsteps / nsubs) +frames = 0:round(Int, nsteps / nsubs) + +record(fig, "barotropicqgql_betaforced.mp4", frames, framerate = 18) do j if j % (1000 / nsubs) == 0 cfl = clock.dt * maximum([maximum(vars.u .+ vars.U) / grid.dx, maximum(vars.v) / grid.dy]) @@ -280,19 +275,23 @@ anim = @animate for j = 0:round(Int, nsteps / nsubs) println(log) end - p[1][1][:z] = Array(@. vars.zeta + vars.Zeta) - p[1][:title] = "vorticity, μt=" * @sprintf("%.2f", μ * clock.t) - p[4][1][:z] = Array(@. vars.psi + vars.Psi) - p[2][1][:x] = Array(mean(vars.Zeta, dims=1)') - p[5][1][:x] = Array(mean(vars.U, dims=1)') - push!(p[3][1], μ * E.t[E.i], E.data[E.i]) - push!(p[6][1], μ * Z.t[Z.i], Z.data[Z.i]) + ζ[] = @. ζ̄ + ζ′ + ψ[] = @. ψ̄ + ψ′ + ζ̄ₘ[] = vec(mean(ζ̄, dims=1)) + ūₘ[] = vec(mean(prob.vars.U, dims=1)) + + μt.val = μ * E.t[1:E.i] + energy[] = E.data[1:E.i] + enstrophy[] = Z.data[1:E.i] + title_ζ[] = @sprintf("vorticity, μt = %.2f", μ * clock.t) + stepforward!(prob, diags, nsubs) BarotropicQGQL.updatevars!(prob) end +nothing # hide -mp4(anim, "barotropicqgql_betaforced.mp4", fps=18) +# ![](barotropicqgql_betaforced.mp4) # ## Save diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 644c8921..68caa332 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -1,4 +1,4 @@ -# # Phillips model of Baroclinic Instability +# # [Phillips model of Baroclinic Instability](@ref multilayerqg_2layer_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/multilayerqg_2layer.ipynb). # @@ -12,13 +12,13 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Plots, Printf" +# pkg"add GeophysicalFlows, CairoMakie, Printf" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, Plots, Printf +using GeophysicalFlows, CairoMakie, Printf using Random: seed! @@ -46,12 +46,12 @@ L = 2π # domain size nlayers = 2 # number of layers f₀, g = 1, 1 # Coriolis parameter and gravitational constant - H = [0.2, 0.8] # the rest depths of each layer - ρ = [4.0, 5.0] # the density of each layer +H = [0.2, 0.8] # the rest depths of each layer +ρ = [4.0, 5.0] # the density of each layer - U = zeros(nlayers) # the imposed mean zonal flow in each layer - U[1] = 1.0 - U[2] = 0.0 +U = zeros(nlayers) # the imposed mean zonal flow in each layer +U[1] = 1.0 +U[2] = 0.0 nothing # hide @@ -61,9 +61,8 @@ nothing # hide # at every time-step that removes enstrophy at high wavenumbers and, thereby, # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. -prob = MultiLayerQG.Problem(nlayers, dev; - nx=n, Lx=L, f₀=f₀, g=g, H=H, ρ=ρ, U=U, μ=μ, β=β, - dt=dt, stepper=stepper, aliased_fraction=0) +prob = MultiLayerQG.Problem(nlayers, dev; nx=n, Lx=L, f₀, g, H, ρ, U, μ, β, + dt, stepper, aliased_fraction=0) nothing # hide # and define some shortcuts. @@ -92,7 +91,7 @@ nothing # hide # ## Diagnostics # Create Diagnostics -- `energies` function is imported at the top. -E = Diagnostic(MultiLayerQG.energies, prob; nsteps=nsteps) +E = Diagnostic(MultiLayerQG.energies, prob; nsteps) diags = [E] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -112,7 +111,9 @@ if !isdir(plotpath); mkdir(plotpath); end nothing # hide # And then create Output -get_sol(prob) = sol # extracts the Fourier-transformed solution + +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution + function get_u(prob) sol, params, vars, grid = prob.sol, prob.params, prob.vars, prob.grid @@ -130,107 +131,136 @@ nothing # hide # ## Visualizing the simulation -# We define a function that plots the potential vorticity field and the evolution of energy -# and enstrophy. Note that when plotting, we decorate the variable to be plotted with `Array()` +# We create a figure using Makie's [`Observable`](https://makie.juliaplots.org/stable/documentation/nodes/)s +# that plots the potential vorticity field and the evolution of energy and enstrophy. +# Note that when plotting, we decorate the variable to be plotted with `Array()` # to make sure it is brought back on the CPU when `vars` live on the GPU. -symlims(data) = maximum(abs.(extrema(data))) |> q -> (-q, q) +Lx, Ly = grid.Lx, grid.Ly -function plot_output(prob) - Lx, Ly = prob.grid.Lx, prob.grid.Ly - - l = @layout Plots.grid(2, 3) - p = plot(layout=l, size = (1000, 600)) - - for m in 1:nlayers - heatmap!(p[(m-1) * 3 + 1], x, y, Array(vars.q[:, :, m]'), - aspectratio = 1, - legend = false, - c = :balance, - xlims = (-Lx/2, Lx/2), - ylims = (-Ly/2, Ly/2), - clims = symlims, - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "q_"*string(m), - framestyle = :box) - - contourf!(p[(m-1) * 3 + 2], x, y, Array(vars.ψ[:, :, m]'), - levels = 8, - aspectratio = 1, - legend = false, - c = :viridis, - xlims = (-Lx/2, Lx/2), - ylims = (-Ly/2, Ly/2), - clims = symlims, - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "ψ_"*string(m), - framestyle = :box) - end +title_KE = Observable(@sprintf("μt = %.2f", μ * clock.t)) - plot!(p[3], 2, - label = ["KE₁" "KE₂"], - legend = :bottomright, - linewidth = 2, - alpha = 0.7, - xlims = (-0.1, 2.35), - ylims = (1e-9, 1e0), - yscale = :log10, - yticks = 10.0.^(-9:0), - xlabel = "μt") - - plot!(p[6], 1, - label = "PE", - legend = :bottomright, - linecolor = :red, - linewidth = 2, - alpha = 0.7, - xlims = (-0.1, 2.35), - ylims = (1e-9, 1e0), - yscale = :log10, - yticks = 10.0.^(-9:0), - xlabel = "μt") +q₁ = Observable(Array(vars.q[:, :, 1])) +ψ₁ = Observable(Array(vars.ψ[:, :, 1])) +q₂ = Observable(Array(vars.q[:, :, 2])) +ψ₂ = Observable(Array(vars.ψ[:, :, 2])) +function compute_levels(maxf, nlevels=8) + ## -max(|f|):...:max(|f|) + levelsf = @lift collect(range(-$maxf, stop = $maxf, length=nlevels)) + + ## only positive + levelsf⁺ = @lift collect(range($maxf/(nlevels-1), stop = $maxf, length=Int(nlevels/2))) + + ## only negative + levelsf⁻ = @lift collect(range(-$maxf, stop = -$maxf/(nlevels-1), length=Int(nlevels/2))) + + return levelsf, levelsf⁺, levelsf⁻ end -nothing # hide + +maxψ₁ = Observable(maximum(abs, vars.ψ[:, :, 1])) +maxψ₂ = Observable(maximum(abs, vars.ψ[:, :, 2])) + +levelsψ₁, levelsψ₁⁺, levelsψ₁⁻ = compute_levels(maxψ₁) +levelsψ₂, levelsψ₂⁺, levelsψ₂⁻ = compute_levels(maxψ₂) + +KE₁ = Observable(Point2f[(μ * E.t[1], E.data[1][1][1])]) +KE₂ = Observable(Point2f[(μ * E.t[1], E.data[1][1][2])]) +PE = Observable(Point2f[(μ * E.t[1], E.data[1][2])]) + +fig = Figure(resolution=(1000, 600)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axq₁ = Axis(fig[1, 1]; title = "q₁", axis_kwargs...) + +axψ₁ = Axis(fig[2, 1]; title = "ψ₁", axis_kwargs...) + +axq₂ = Axis(fig[1, 2]; title = "q₂", axis_kwargs...) + +axψ₂ = Axis(fig[2, 2]; title = "ψ₂", axis_kwargs...) + +axKE = Axis(fig[1, 3], + xlabel = "μ t", + ylabel = "KE", + title = title_KE, + yscale = log10, + limits = ((-0.1, 2.6), (1e-9, 5))) + +axPE = Axis(fig[2, 3], + xlabel = "μ t", + ylabel = "PE", + yscale = log10, + limits = ((-0.1, 2.6), (1e-9, 5))) + +heatmap!(axq₁, x, y, q₁; colormap = :balance) + +heatmap!(axq₂, x, y, q₂; colormap = :balance) + +contourf!(axψ₁, x, y, ψ₁; + levels = levelsψ₁, colormap = :viridis, extendlow = :auto, extendhigh = :auto) + contour!(axψ₁, x, y, ψ₁; + levels = levelsψ₁⁺, color=:black) + contour!(axψ₁, x, y, ψ₁; + levels = levelsψ₁⁻, color=:black, linestyle = :dash) + +contourf!(axψ₂, x, y, ψ₂; + levels = levelsψ₂, colormap = :viridis, extendlow = :auto, extendhigh = :auto) + contour!(axψ₂, x, y, ψ₂; + levels = levelsψ₂⁺, color=:black) + contour!(axψ₂, x, y, ψ₂; + levels = levelsψ₂⁻, color=:black, linestyle = :dash) + +ke₁ = lines!(axKE, KE₁; linewidth = 3) +ke₂ = lines!(axKE, KE₂; linewidth = 3) +Legend(fig[1, 4], [ke₁, ke₂,], ["KE₁", "KE₂"]) + +lines!(axPE, PE; linewidth = 3) + +fig # ## Time-stepping the `Problem` forward # Finally, we time-step the `Problem` forward in time. -p = plot_output(prob) - startwalltime = time() -anim = @animate for j = 0:round(Int, nsteps / nsubs) +frames = 0:round(Int, nsteps / nsubs) + +record(fig, "multilayerqg_2layer.mp4", frames, framerate = 18) do j if j % (1000 / nsubs) == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) - log = @sprintf("step: %04d, t: %.1f, cfl: %.2f, KE₁: %.3e, KE₂: %.3e, PE: %.3e, walltime: %.2f min", clock.step, clock.t, cfl, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60) + log = @sprintf("step: %04d, t: %.1f, cfl: %.2f, KE₁: %.3e, KE₂: %.3e, PE: %.3e, walltime: %.2f min", + clock.step, clock.t, cfl, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60) println(log) end - for m in 1:nlayers - p[(m-1) * 3 + 1][1][:z] = Array(vars.q[:, :, m]) - p[(m-1) * 3 + 2][1][:z] = Array(vars.ψ[:, :, m]) - end - - push!(p[3][1], μ * E.t[E.i], E.data[E.i][1][1]) - push!(p[3][2], μ * E.t[E.i], E.data[E.i][1][2]) - push!(p[6][1], μ * E.t[E.i], E.data[E.i][2][1]) + q₁[] = vars.q[:, :, 1] + ψ₁[] = vars.ψ[:, :, 1] + q₂[] = vars.q[:, :, 2] + ψ₂[] = vars.ψ[:, :, 2] + + maxψ₁[] = maximum(abs, vars.ψ[:, :, 1]) + maxψ₂[] = maximum(abs, vars.ψ[:, :, 2]) + + KE₁[] = push!(KE₁[], Point2f(μ * E.t[E.i], E.data[E.i][1][1])) + KE₂[] = push!(KE₂[], Point2f(μ * E.t[E.i], E.data[E.i][1][2])) + PE[] = push!(PE[] , Point2f(μ * E.t[E.i], E.data[E.i][2])) + + title_KE[] = @sprintf("μ t = %.2f", μ * clock.t) stepforward!(prob, diags, nsubs) MultiLayerQG.updatevars!(prob) end +nothing # hide -mp4(anim, "multilayerqg_2layer.mp4", fps=18) +# ![](multilayerqg_2layer.mp4) # ## Save diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index ddbfd7e2..3d275ce0 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -1,4 +1,4 @@ -# # Decaying barotropic QG beta-plane turbulence +# # [Decaying barotropic QG beta-plane turbulence](@id singlelayerqg_betadecay_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_betadecay.ipynb). # @@ -10,13 +10,13 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Plots, Printf, Statistics, Random" +# pkg"add GeophysicalFlows, CairoMakie, Printf, Statistics, Random" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, Plots, Printf, Random +using GeophysicalFlows, CairoMakie, Printf, Random using Statistics: mean @@ -51,8 +51,7 @@ nothing # hide # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. # Thus, we choose not to do any dealiasing by providing `aliased_fraction=0`. -prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, - dt=dt, stepper=stepper, aliased_fraction=0) +prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, β, μ, dt, stepper, aliased_fraction=0) nothing # hide # and define some shortcuts @@ -87,42 +86,35 @@ nothing # hide # the variable to be plotted with `Array()` to make sure it is brought back on the CPU when # `vars` live on the GPU. -p1 = heatmap(x, y, Array(vars.q'), - aspectratio = 1, - c = :balance, - clim = (-12, 12), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "initial vorticity ∂v/∂x-∂u/∂y", - framestyle = :box) - -p2 = contourf(x, y, Array(vars.ψ'), - aspectratio = 1, - c = :viridis, - levels = range(-0.7, stop=0.7, length=20), - clim = (-0.35, 0.35), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "initial streamfunction ψ", - framestyle = :box) - -layout = @layout Plots.grid(1, 2) -p = plot(p1, p2, layout = layout, size = (800, 360)) +fig = Figure(resolution = (800, 360)) +axq = Axis(fig[1, 1]; + xlabel = "x", + ylabel = "y", + title = "initial vorticity ∂v/∂x-∂u/∂y", + aspect = 1, + limits = ((-grid.Lx/2, grid.Lx/2), (-grid.Ly/2, grid.Ly/2)) + ) + +axψ = Axis(fig[1, 2]; + xlabel = "x", + ylabel = "y", + title = "initial streamfunction ψ", + aspect = 1, + limits = ((-grid.Lx/2, grid.Lx/2), (-grid.Ly/2, grid.Ly/2)) + ) + +heatmap!(axq, x, y, Array(vars.q); colormap = :balance) + +contourf!(axψ, x, y, Array(vars.ψ); colormap = :viridis) + +fig # ## Diagnostics # Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top. -E = Diagnostic(SingleLayerQG.energy, prob; nsteps=nsteps) -Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps=nsteps) +E = Diagnostic(SingleLayerQG.energy, prob; nsteps) +Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps) diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -142,89 +134,76 @@ if !isdir(plotpath); mkdir(plotpath); end nothing # hide # and then create Output. -get_sol(prob) = sol # extracts the Fourier-transformed solution +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution out = Output(prob, filename, (:sol, get_sol)) nothing # hide # ## Visualizing the simulation -# We define a function that plots the vorticity and streamfunction and -# their corresponding zonal mean structure. - -function plot_output(prob) - q = Array(prob.vars.q) - ψ = Array(prob.vars.ψ) - q̄ = Array(mean(q, dims=1)') - ū = Array(mean(prob.vars.u, dims=1)') - - pq = heatmap(x, y, q', - aspectratio = 1, - legend = false, - c = :balance, - clim = (-12, 12), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "vorticity ∂v/∂x-∂u/∂y", - framestyle = :box) - - pψ = contourf(x, y, ψ', - aspectratio = 1, - legend = false, - c = :viridis, - levels = range(-0.7, stop=0.7, length=20), - clim = (-0.35, 0.35), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "streamfunction ψ", - framestyle = :box) - - pqm = plot(q̄, y, - legend = false, - linewidth = 2, - alpha = 0.7, - yticks = -3:3, - xlims = (-2.2, 2.2), - xlabel = "zonal mean q", - ylabel = "y") - plot!(pqm, 0*y, y, linestyle=:dash, linecolor=:black) - - pum = plot(ū, y, - legend = false, - linewidth = 2, - alpha = 0.7, - yticks = -3:3, - xlims = (-0.55, 0.55), - xlabel = "zonal mean u", - ylabel = "y") - plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black) - - layout = @layout Plots.grid(2, 2) - p = plot(pq, pqm, pψ, pum, layout = layout, size = (800, 720)) - - return p -end -nothing # hide +# We plot the vorticity and streamfunction and their corresponding zonal mean structure. + +Lx, Ly = grid.Lx, grid.Ly + +title_q = Observable(@sprintf("vorticity, t = %.2f", clock.t)) +title_ψ = "streamfunction ψ" + +fig = Figure(resolution=(800, 720)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axq = Axis(fig[1, 1]; title = title_q, axis_kwargs...) + +axψ = Axis(fig[2, 1]; title = title_ψ, axis_kwargs...) + +axq̄ = Axis(fig[1, 2], + xlabel = "zonal mean vorticity", + ylabel = "y", + aspect = 1, + limits = ((-2.1, 2.1), (-Ly/2, Ly/2))) + +axū = Axis(fig[2, 2], + xlabel = "zonal mean u", + ylabel = "y", + aspect = 1, + limits = ((-0.5, 0.5), (-Ly/2, Ly/2))) + +q = Observable(Array(vars.q)) +ψ = Observable(Array(vars.ψ)) +q̄ₘ = Observable(Array(vec(mean(vars.q, dims=1)))) +ūₘ = Observable(Array(vec(mean(vars.u, dims=1)))) + +heatmap!(axq, x, y, q; + colormap = :balance, colorrange = (-12, 12)) + +levels = collect(range(-0.7, stop=0.7, length=20)) + +contourf!(axψ, x, y, ψ; + levels, colormap = :viridis, colorrange = (-0.35, 0.35)) +contour!(axψ, x, y, ψ; + levels, color = :black) + +lines!(axq̄, q̄ₘ, y; linewidth = 3) +lines!(axq̄, 0y, y; linewidth = 1, linestyle = :dash) + +lines!(axū, ūₘ, y; linewidth = 3) +lines!(axū, 0y, y; lindewidth = 1, linestyle = :dash) + +fig # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() -p = plot_output(prob) - -anim = @animate for j = 0:round(Int, nsteps/nsubs) +frames = 0:round(Int, nsteps / nsubs) +record(fig, "singlelayerqg_betadecay.mp4", frames, framerate = 8) do j if j % round(Int, nsteps/nsubs / 4) == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) @@ -232,20 +211,22 @@ anim = @animate for j = 0:round(Int, nsteps/nsubs) clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60) println(log) - end + end + + q[] = vars.q + ψ[] = vars.ψ + q̄ₘ[] = vec(mean(vars.q, dims=1)) + ūₘ[] = vec(mean(vars.u, dims=1)) - p[1][1][:z] = Array(vars.q) - p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t) - p[3][1][:z] = Array(vars.ψ) - p[2][1][:x] = Array(mean(vars.q, dims=1)') - p[4][1][:x] = Array(mean(vars.u, dims=1)') + title_q[] = @sprintf("vorticity, t = %.2f", clock.t) stepforward!(prob, diags, nsubs) SingleLayerQG.updatevars!(prob) - end +nothing # hide + +# ![](singlelayerqg_betadecay.mp4) -mp4(anim, "singlelayerqg_betadecay.mp4", fps=8) # ## Save diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 7e19e26e..bbceac9d 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -1,4 +1,4 @@ -# # Forced-dissipative barotropic QG beta-plane turbulence +# # [Forced-dissipative barotropic QG beta-plane turbulence](@id singlelayerqg_betaforced_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_betaforced.ipynb). # @@ -11,16 +11,19 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, CUDA, Random, Printf, Plots, Statistics" +# pkg"add GeophysicalFlows, CUDA, JLD2, CairoMakie, Statistics" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, CUDA, Random, Printf, Plots +using GeophysicalFlows, CUDA, JLD2, CairoMakie, Random, Printf using Statistics: mean +using LinearAlgebra: ldiv! + parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU @@ -33,8 +36,9 @@ nothing # hide n = 128 # 2D resolution: n² grid points stepper = "FilteredRK4" # timestepper dt = 0.05 # timestep - nsteps = 8000 # total number of time-steps - nsubs = 10 # number of time-steps for intermediate logging/plotting (nsteps must be multiple of nsubs) + nsteps = 8000 # total number of timesteps + save_substeps = 10 # number of timesteps after which output is saved + nothing # hide @@ -94,13 +98,14 @@ nothing # hide # We use `stepper = "FilteredRK4"`. Filtered timesteppers apply a wavenumber-filter # at every time-step that removes enstrophy at high wavenumbers and, thereby, # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. -prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, β=β, μ=μ, dt=dt, stepper=stepper, +prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, β, μ, dt, stepper, calcF=calcF!, stochastic=true) nothing # hide # Let's define some shortcuts. sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -x, y = grid.x, grid.y +x, y = grid.x, grid.y +Lx, Ly = grid.Lx, grid.Ly nothing # hide @@ -109,19 +114,20 @@ nothing # hide # `vars` live on the GPU. calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) -heatmap(x, y, Array(irfft(vars.Fh, grid.nx)'), - aspectratio = 1, - c = :balance, - clim = (-8, 8), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, +fig = Figure() + +ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", - title = "a forcing realization", - framestyle = :box) - + aspect = 1, + title = "a forcing realization", + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); + colormap = :balance, colorrange = (-8, 8)) + +fig + # ## Setting initial conditions @@ -132,8 +138,8 @@ SingleLayerQG.set_q!(prob, ArrayType(dev)(zeros(grid.nx, grid.ny))) # ## Diagnostics # Create Diagnostic -- `energy` and `enstrophy` are functions imported at the top. -E = Diagnostic(SingleLayerQG.energy, prob; nsteps=nsteps) -Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps=nsteps) +E = Diagnostic(SingleLayerQG.energy, prob; nsteps, freq=save_substeps) +Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps, freq=save_substeps) diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -144,7 +150,7 @@ nothing # hide filepath = "." plotpath = "./plots_forcedbetaturb" plotname = "snapshots" -filename = joinpath(filepath, "forcedbetaturb.jld2") +filename = joinpath(filepath, "singlelayerqg_forcedbeta.jld2") nothing # hide # Do some basic file management, @@ -153,99 +159,28 @@ if !isdir(plotpath); mkdir(plotpath); end nothing # hide # and then create Output. -get_sol(prob) = sol # extracts the Fourier-transformed solution -get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx) -out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) -nothing # hide +get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution +function get_u(prob) + vars, grid, sol = prob.vars, prob.grid, prob.sol + + @. vars.qh = sol -# ## Visualizing the simulation + SingleLayerQG.streamfunctionfrompv!(vars.ψh, vars.qh, params, grid) -# We define a function that plots the vorticity and streamfunction fields, their -# corresponding zonal mean structure and timeseries of energy and enstrophy. + ldiv!(vars.u, grid.rfftplan, -im * grid.l .* vars.ψh) -function plot_output(prob) - q = prob.vars.q - ψ = prob.vars.ψ - q̄ = mean(q, dims=1)' - ū = mean(prob.vars.u, dims=1)' - - pq = heatmap(x, y, Array(q'), - aspectratio = 1, - legend = false, - c = :balance, - clim = (-8, 8), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "vorticity ∂v/∂x-∂u/∂y", - framestyle = :box) - - pψ = contourf(x, y, Array(ψ'), - levels = -0.32:0.04:0.32, - aspectratio = 1, - linewidth = 1, - legend = false, - clim = (-0.22, 0.22), - c = :viridis, - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "streamfunction ψ", - framestyle = :box) - - pqm = plot(Array(q̄), y, - legend = false, - linewidth = 2, - alpha = 0.7, - yticks = -3:3, - xlims = (-3, 3), - xlabel = "zonal mean q", - ylabel = "y") - plot!(pqm, 0*y, y, linestyle=:dash, linecolor=:black) - - pum = plot(Array(ū), y, - legend = false, - linewidth = 2, - alpha = 0.7, - yticks = -3:3, - xlims = (-0.5, 0.5), - xlabel = "zonal mean u", - ylabel = "y") - plot!(pum, 0*y, y, linestyle=:dash, linecolor=:black) - - pE = plot(1, - label = "energy", - legend = :bottomright, - linewidth = 2, - alpha = 0.7, - xlims = (-0.1, 4.1), - ylims = (0, 0.05), - xlabel = "μt") - - pZ = plot(1, - label = "enstrophy", - linecolor = :red, - legend = :bottomright, - linewidth = 2, - alpha = 0.7, - xlims = (-0.1, 4.1), - ylims = (0, 3), - xlabel = "μt") - - l = @layout Plots.grid(2, 3) - p = plot(pq, pqm, pE, pψ, pum, pZ, layout=l, size = (1000, 600)) - - return p + return Array(vars.u) end + +output = Output(prob, filename, (:qh, get_sol), (:u, get_u)) nothing # hide +# We first save the problem's grid and other parameters so we can use them later. +saveproblem(output) + +# and then call `saveoutput(output)` once to save the initial state. +saveoutput(output) # ## Time-stepping the `Problem` forward @@ -253,38 +188,130 @@ nothing # hide startwalltime = time() -p = plot_output(prob) - -anim = @animate for j = 0:Int(nsteps / nsubs) - - if j % (1000 / nsubs) == 0 +while clock.step <= nsteps + if clock.step % 50save_substeps == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, walltime: %.2f min", clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60) println(log) - end - - p[1][1][:z] = Array(vars.q) - p[1][:title] = "vorticity, μt="*@sprintf("%.2f", μ * clock.t) - p[4][1][:z] = Array(vars.ψ) - p[2][1][:x] = Array(mean(vars.q, dims=1)') - p[5][1][:x] = Array(mean(vars.u, dims=1)') - push!(p[3][1], μ * E.t[E.i], E.data[E.i]) - push!(p[6][1], μ * Z.t[Z.i], Z.data[Z.i]) - - stepforward!(prob, diags, nsubs) + end + + stepforward!(prob, diags, save_substeps) SingleLayerQG.updatevars!(prob) + + if clock.step % save_substeps == 0 + saveoutput(output) + end end -mp4(anim, "singlelayerqg_betaforced.mp4", fps=18) +savediagnostic(E, "energy", output.path) +savediagnostic(Z, "enstrophy", output.path) -# ## Save +# ## Load saved output and visualize -# Finally, we can save, e.g., the last snapshot via -# ```julia -# savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step) -# savefig(savename) -# ``` +# We now have output from our simulation saved in `singlelayerqg_forcedbeta.jld2` which +# we can load to create a time series for the fields we are interested in. + +file = jldopen(output.path) + +iterations = parse.(Int, keys(file["snapshots/t"])) +t = [file["snapshots/t/$i"] for i ∈ iterations] + +qh = [file["snapshots/qh/$i"] for i ∈ iterations] +u = [file["snapshots/u/$i"] for i ∈ iterations] + +E_t = file["diagnostics/energy/t"] +Z_t = file["diagnostics/enstrophy/t"] +E_data = file["diagnostics/energy/data"] +Z_data = file["diagnostics/enstrophy/data"] + +x, y = file["grid/x"], file["grid/y"] +nx, ny = file["grid/nx"], file["grid/ny"] +Lx, Ly = file["grid/Lx"], file["grid/Ly"] + +# We create a figure using Makie's [`Observable`](https://makie.juliaplots.org/stable/documentation/nodes/)s + +j = Observable(1) + +q = @lift irfft(qh[$j], nx) +ψ = @lift irfft(- Array(grid.invKrsq) .* qh[$j], nx) +q̄ = @lift real(ifft(qh[$j][1, :] / ny)) +ū = @lift vec(mean(u[$j], dims=1)) + +title_q = @lift @sprintf("vorticity, μt = %.2f", μ * t[$j]) + +energy = Observable([Point2f(E_t[1], E_data[1])]) +enstrophy = Observable([Point2f(Z_t[1], Z_data[1])]) + +fig = Figure(resolution=(1000, 600)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axq = Axis(fig[1, 1]; title = title_q, axis_kwargs...) + +axψ = Axis(fig[2, 1]; title = "streamfunction ψ", axis_kwargs...) + +axq̄ = Axis(fig[1, 2], + xlabel = "zonal mean vorticity", + ylabel = "y", + aspect = 1, + limits = ((-3, 3), (-Ly/2, Ly/2))) + +axū = Axis(fig[2, 2], + xlabel = "zonal mean u", + ylabel = "y", + aspect = 1, + limits = ((-0.5, 0.5), (-Ly/2, Ly/2))) + +axE = Axis(fig[1, 3], + xlabel = "μ t", + ylabel = "energy", + aspect = 1, + limits = ((-0.1, 4.1), (0, 0.055))) + +axZ = Axis(fig[2, 3], + xlabel = "μ t", + ylabel = "enstrophy", + aspect = 1, + limits = ((-0.1, 4.1), (0, 3.1))) + +heatmap!(axq, x, y, q; + colormap = :balance, colorrange = (-8, 8)) + +levels = collect(-0.32:0.04:0.32) + +contourf!(axψ, x, y, ψ; + levels, colormap = :viridis, colorrange = (-0.22, 0.22)) +contour!(axψ, x, y, ψ; + levels, color = :black) + +lines!(axq̄, q̄, y; linewidth = 3) +lines!(axq̄, 0y, y; linewidth = 1, linestyle=:dash) + +lines!(axū, ū, y; linewidth = 3) +lines!(axū, 0y, y; linewidth = 1, linestyle=:dash) + +lines!(axE, energy; linewidth = 3) +lines!(axZ, enstrophy; linewidth = 3, color = :red) + +fig + + +# We are now ready to animate all saved output. + +frames = 1:length(t) +record(fig, "singlelayerqg_betaforced.mp4", frames, framerate = 18) do i + j[] = i + + energy[] = push!(energy[], Point2f(μ * E_t[i], E_data[i])) + enstrophy[] = push!(enstrophy[], Point2f(μ * Z_t[i], Z_data[i])) +end +nothing # hide + +# ![](singlelayerqg_betaforced.mp4) diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index b2e01ca2..b2775673 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -1,4 +1,4 @@ -# # SingleLayerQG decaying 2D turbulence with and without finite Rossby radius of deformation +# # [SingleLayerQG decaying 2D turbulence with and without finite Rossby radius of deformation](@id singlelayerqg_decaying_barotropic_equivalentbarotropic_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_decaying_barotropic_equivalentbarotropic.ipynb). # @@ -11,15 +11,16 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Printf, Random, Plots" +# pkg"add GeophysicalFlows, Printf, Random, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, Printf, Random, Plots +using GeophysicalFlows, Printf, Random, CairoMakie using GeophysicalFlows: peakedisotropicspectrum +using LinearAlgebra: ldiv! using Random: seed! @@ -54,10 +55,10 @@ nothing # hide # thereby, stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. # Thus, we choose not to do any dealiasing by providing `aliased_fraction=0`. -prob_bqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, - dt=dt, stepper="FilteredRK4", aliased_fraction=0) -prob_eqbqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, deformation_radius = deformation_radius, - dt=dt, stepper="FilteredRK4", aliased_fraction=0) +stepper="FilteredRK4" + +prob_bqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, dt, stepper, aliased_fraction=0) +prob_eqbqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, deformation_radius, dt, stepper, aliased_fraction=0) nothing # hide @@ -91,43 +92,47 @@ SingleLayerQG.set_q!(prob_eqbqg, q₀_eqbqg) nothing # hide -# Let's plot the initial vorticity field for each problem. A function that returns relative -# vorticity from each problem's state variable will prove useful. Note that when plotting, we -# decorate the variable to be plotted with `Array()` to make sure it is brought back on the -# CPU when the variable lives on the GPU. -relativevorticity(prob) = irfft(-prob.grid.Krsq .* prob.vars.ψh, prob.grid.nx) - -x, y = prob_bqg.grid.x, prob_bqg.grid.y - -p_bqg = heatmap(x, y, Array(relativevorticity(prob_bqg)'), - 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 = "barotropic\n ∇²ψ, t=" * @sprintf("%.2f", prob_bqg.clock.t), - framestyle = :box) - -p_eqbqg = heatmap(x, y, Array(relativevorticity(prob_eqbqg)'), - 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 = "equivalent barotropic; deformation radius: " * @sprintf("%.2f", prob_eqbqg.params.deformation_radius) * "\n ∇²ψ, t=" * @sprintf("%.2f", prob_eqbqg.clock.t), - framestyle = :box) - -l = @layout Plots.grid(1, 2) -p = plot(p_bqg, p_eqbqg, layout = l, size = (800, 380)) +# Let's plot the initial vorticity field for each problem. Note that when plotting, we decorate +# the variable to be plotted with `Array()` to make sure it is brought back on the CPU when +# `vars` live on the GPU. + +function relativevorticity(prob) + vars, grid = prob.vars, prob.grid + + ldiv!(vars.q, grid.rfftplan, - grid.Krsq .* vars.ψh) + + return vars.q +end + +x, y = prob_bqg.grid.x, prob_bqg.grid.y +Lx, Ly = prob_bqg.grid.Lx, prob_bqg.grid.Ly + +fig = Figure(resolution=(800, 380)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +t_bqg = Observable(prob_bqg.clock.t) +t_eqbqg = Observable(prob_eqbqg.clock.t) +title_bqg = @lift "barotropic\n ∇²ψ, t=" * @sprintf("%.2f", $t_bqg) +title_eqbqg = @lift "equivalent barotropic; deformation radius: " * @sprintf("%.2f", prob_eqbqg.params.deformation_radius) * "\n ∇²ψ, t=" * @sprintf("%.2f", $t_eqbqg) + +ax1 = Axis(fig[1, 1]; title = title_bqg, axis_kwargs...) +ax2 = Axis(fig[1, 2]; title = title_eqbqg, axis_kwargs...) + +ζ_bqg = Observable(Array(relativevorticity(prob_bqg))) +ζ_eqbqg = Observable(Array(relativevorticity(prob_eqbqg))) + +heatmap!(ax1, x, y, ζ_bqg; + colormap = :balance, colorrange = (-40, 40)) + +heatmap!(ax2, x, y, ζ_eqbqg; + colormap = :balance, colorrange = (-40, 40)) + +fig # ## Time-stepping the `Problem` forward @@ -137,7 +142,7 @@ startwalltime = time() cfl(prob) = prob.clock.dt * maximum([maximum(prob.vars.u) / prob.grid.dx, maximum(prob.vars.v) / prob.grid.dy]) -anim = @animate for j = 0:Int(nsteps/nsubs) +record(fig, "singlelayerqg_barotropic_equivalentbarotropic.mp4", 0:Int(nsteps/nsubs), framerate = 18) do j if j % (1000 / nsubs) == 0 log_bqg = @sprintf("barotropic; step: %04d, t: %d, cfl: %.2f, walltime: %.2f min", prob_bqg.clock.step, prob_bqg.clock.t, cfl(prob_bqg), (time()-startwalltime)/60) @@ -148,16 +153,17 @@ anim = @animate for j = 0:Int(nsteps/nsubs) println(log_eqbqg) end - p[1][1][:z] = Array(relativevorticity(prob_bqg)) - p[1][:title] = "barotropic\n ∇²ψ, t=" * @sprintf("%.2f", prob_bqg.clock.t) - p[2][1][:z] = Array(relativevorticity(prob_eqbqg)) - p[2][:title] = "equivalent barotropic; deformation radius: " * @sprintf("%.2f", prob_eqbqg.params.deformation_radius) * "\n ∇²ψ, t=" * @sprintf("%.2f", prob_eqbqg.clock.t) - stepforward!(prob_bqg, nsubs) SingleLayerQG.updatevars!(prob_bqg) stepforward!(prob_eqbqg, nsubs) SingleLayerQG.updatevars!(prob_eqbqg) + + t_bqg[] = prob_bqg.clock.t + t_eqbqg[] = prob_eqbqg.clock.t + ζ_bqg[] = relativevorticity(prob_bqg) + ζ_eqbqg[] = relativevorticity(prob_eqbqg) end +nothing # hide -mp4(anim, "singlelayerqg_barotropic_equivalentbarotropic.mp4", fps=18) +# ![](singlelayerqg_barotropic_equivalentbarotropic.mp4) diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 1de9b047..0a320ba8 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -1,4 +1,4 @@ -# # Decaying barotropic QG turbulence over topography +# # [Decaying barotropic QG turbulence over topography](@id singlelayerqg_decaying_topography) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/singlelayerqg_decaying_topography.ipynb). # @@ -10,13 +10,13 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Plots, Printf, Random, Statistics" +# pkg"add GeophysicalFlows, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, Plots, Printf, Random +using GeophysicalFlows, CairoMakie, Printf, Random using Statistics: mean @@ -45,7 +45,7 @@ nothing # hide # Define the topographic potential vorticity, ``\eta = f_0 h(x, y)/H``. The topography here is # an elliptical mount at ``(x, y) = (1, 1)``, and an elliptical depression at ``(x, y) = (-1, -1)``. σx, σy = 0.4, 0.8 -topographicPV(x, y) = 3exp(-(x-1)^2/(2σx^2) -(y-1)^2/(2σy^2)) - 2exp(-(x+1)^2/(2σx^2) -(y+1)^2/(2σy^2)) +topographicPV(x, y) = 3exp(-(x - 1)^2 / 2σx^2 - (y - 1)^2 / 2σy^2) - 2exp(- (x + 1)^2 / 2σx^2 - (y + 1)^2 / 2σy^2) nothing # hide # ## Problem setup @@ -57,31 +57,32 @@ nothing # hide # # The topophic PV is prescribed via keyword argument `eta`. prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, eta=topographicPV, - dt=dt, stepper=stepper, aliased_fraction=0) + dt, stepper, aliased_fraction=0) nothing # hide # and define some shortcuts sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -x, y = grid.x, grid.y +x, y = grid.x, grid.y +Lx, Ly = grid.Lx, grid.Ly nothing # hide # and let's plot the topographic PV. Note that when plotting, we decorate the variable to be # plotted with `Array()` to make sure it is brought back on the CPU when the variable lives # on the GPU. -contourf(grid.x, grid.y, Array(params.eta'), - aspectratio = 1, - linewidth = 0, - levels = 10, - c = :balance, - clim = (-3, 3), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "topographic PV η=f₀h/H") +η = Array(params.eta) + +fig = Figure() +ax = Axis(fig[1, 1]; + xlabel = "x", + ylabel = "y", + title = "topographic PV η=f₀h/H", + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +contourf!(ax, x, y, η; + levels = collect(-3:0.4:3), colormap = :balance, colorrange = (-3, 3)) + +fig # ## Setting initial conditions @@ -104,44 +105,47 @@ qi = irfft(qih, grid.nx) SingleLayerQG.set_q!(prob, qi) nothing # hide -# Let's plot the initial vorticity field: - -p1 = heatmap(x, y, Array(vars.q'), - aspectratio = 1, - c = :balance, - clim = (-8, 8), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "initial vorticity ∂v/∂x-∂u/∂y", - framestyle = :box) - -p2 = contourf(x, y, Array(vars.ψ'), - aspectratio = 1, - c = :viridis, - levels = range(-0.25, stop=0.25, length=11), - clim = (-0.25, 0.25), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "initial streamfunction ψ", - framestyle = :box) - -layout = @layout Plots.grid(1, 2) -p = plot(p1, p2, layout=layout, size = (800, 360)) +# Let's plot the initial vorticity and streamfunction. + +q = Observable(Array(vars.q)) +ψ = Observable(Array(vars.ψ)) + +fig = Figure(resolution=(800, 380)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +title_q = Observable("initial vorticity ∂v/∂x-∂u/∂y") +axq = Axis(fig[1, 1]; title = title_q, axis_kwargs...) + +title_ψ = Observable("initial streamfunction ψ") +axψ = Axis(fig[1, 3]; title = title_ψ, axis_kwargs...) + +hm = heatmap!(axq, x, y, q; + colormap = :balance, colorrange = (-8, 8)) + +Colorbar(fig[1, 2], hm) + +levels = collect(range(-0.28, stop=0.28, length=11)) + +hc = contourf!(axψ, x, y, ψ; + levels, colormap = :viridis, colorrange = (-0.28, 0.28), + extendlow = :auto, extendhigh = :auto) +contour!(axψ, x, y, ψ; + levels, color = :black) + +Colorbar(fig[1, 4], hc) + +fig # ## Diagnostics # Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top. -E = Diagnostic(SingleLayerQG.energy, prob; nsteps=nsteps) -Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps=nsteps) +E = Diagnostic(SingleLayerQG.energy, prob; nsteps) +Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps) diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -158,76 +162,35 @@ if isfile(filename); rm(filename); end nothing # hide # and then create Output. -get_sol(prob) = sol # extracts the Fourier-transformed solution +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution out = Output(prob, filename, (:sol, get_sol)) nothing # hide # ## Visualizing the simulation -# We define a function that plots the vorticity and streamfunction and -# their corresponding zonal mean structure. - -function plot_output(prob) - q = prob.vars.q - ψ = prob.vars.ψ - η = prob.params.eta - - pq = heatmap(x, y, Array(q'), - aspectratio = 1, - legend = false, - c = :balance, - clim = (-6, 6), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "vorticity ∂v/∂x-∂u/∂y", - framestyle = :box) - - contour!(pq, x, y, Array(η'), - levels=0.5:0.5:3, - lw=2, c=:black, ls=:solid, alpha=0.7) - - contour!(pq, x, y, Array(η'), - levels=-2:0.5:-0.5, - lw=2, c=:black, ls=:dash, alpha=0.7) - - pψ = contourf(x, y, Array(ψ'), - aspectratio = 1, - legend = false, - c = :viridis, - levels = range(-0.75, stop=0.75, length=31), - clim = (-0.75, 0.75), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "streamfunction ψ", - framestyle = :box) - - l = @layout Plots.grid(1, 2) - p = plot(pq, pψ, layout = l, size = (800, 360)) - - return p -end +# We modify the figure with the initial state slightly by adding the topography contours +# and mark the time in the title. + +contour!(axq, x, y, η; + levels = collect(0.5:0.5:3), linewidth = 2, color = (:black, 0.5)) + +contour!(axq, x, y, η; + levels = collect(-2:0.5:-0.5), linewidth = 2, color = (:grey, 0.7), linestyle = :dash) + +title_q[] = "vorticity, t=" * @sprintf("%.2f", clock.t) +title_ψ[] = "streamfunction ψ" + nothing # hide # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() -p = plot_output(prob) - -anim = @animate for j = 0:round(Int, nsteps/nsubs) - +record(fig, "singlelayerqg_decaying_topography.mp4", 0:round(Int, nsteps/nsubs), framerate = 12) do j if j % (1000 / nsubs) == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) @@ -235,14 +198,17 @@ anim = @animate for j = 0:round(Int, nsteps/nsubs) clock.step, clock.t, cfl, E.data[E.i], Z.data[Z.i], (time()-startwalltime)/60) println(log) - end + end + + q[] = vars.q + ψ[] = vars.ψ - p[1][1][:z] = Array(vars.q) - p[1][:title] = "vorticity, t="*@sprintf("%.2f", clock.t) - p[2][1][:z] = Array(vars.ψ) + title_q[] = "vorticity, t="*@sprintf("%.2f", clock.t) + title_ψ[] = "streamfunction ψ" stepforward!(prob, diags, nsubs) SingleLayerQG.updatevars!(prob) end +nothing # hide -mp4(anim, "singlelayerqg_decaying_topography.mp4", fps=12) +# ![](singlelayerqg_decaying_topography.mp4) diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index 492f140a..5e95accb 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -1,4 +1,4 @@ -# # Decaying Surface QG turbulence +# # [Decaying Surface QG turbulence](@id surfaceqg_decaying_example) # #md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/literated/surfaceqg_decaying.ipynb). #md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/surfaceqg_decaying.ipynb). @@ -15,13 +15,13 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Plots, Printf, Random, Statistics" +# pkg"add GeophysicalFlows, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, Plots, Printf, Random +using GeophysicalFlows, CairoMakie, Printf, Random using Statistics: mean using Random: seed! @@ -58,12 +58,13 @@ nothing # hide # at every time-step that removes enstrophy at high wavenumbers and, thereby, # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. -prob = SurfaceQG.Problem(dev; nx=n, Lx=L, dt=dt, stepper=stepper, ν=ν, nν=nν) +prob = SurfaceQG.Problem(dev; nx=n, Lx=L, dt, stepper, ν, nν) nothing # hide # Let's define some shortcuts. sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -x, y = grid.x, grid.y +x, y = grid.x, grid.y +Lx, Ly = grid.Lx, grid.Ly #md nothing # hide @@ -71,34 +72,38 @@ x, y = grid.x, grid.y # # We initialize the buoyancy equation with an elliptical vortex. X, Y = gridpoints(grid) -b₀ = @. exp(-(X^2 + 4*Y^2)) +b₀ = @. exp(-(X^2 + 4Y^2)) SurfaceQG.set_b!(prob, b₀) nothing # hide # Let's plot the initial condition. Note that when plotting, we decorate the variable to be # plotted with `Array()` to make sure it is brought back on the CPU when `vars` live on the GPU. -heatmap(x, y, Array(vars.b'), - aspectratio = 1, - c = :deep, - clim = (0, 1), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, + +fig = Figure(resolution = (500, 500)) + +ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", - title = "buoyancy bₛ", - framestyle = :box) + aspect = 1, + title = "buoyancy bₛ", + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +hm = heatmap!(ax, x, y, Array(vars.b); + colormap = :deep, colorrange = (0, 1)) + +Colorbar(fig[1, 2], hm) + +fig # ## Diagnostics # Create Diagnostics; `buoyancy_variance`, `kinetic_energy` and `buoyancy_dissipation` # functions were imported at the top. -B = Diagnostic(SurfaceQG.buoyancy_variance, prob; nsteps=nsteps) -KE = Diagnostic(SurfaceQG.kinetic_energy, prob; nsteps=nsteps) -Dᵇ = Diagnostic(SurfaceQG.buoyancy_dissipation, prob; nsteps=nsteps) +B = Diagnostic(SurfaceQG.buoyancy_variance, prob; nsteps) +KE = Diagnostic(SurfaceQG.kinetic_energy, prob; nsteps) +Dᵇ = Diagnostic(SurfaceQG.buoyancy_dissipation, prob; nsteps) diags = [B, KE, Dᵇ] # A list of Diagnostics types passed to `stepforward!`. Diagnostics are updated every timestep. nothing # hidenothing # hide @@ -122,8 +127,9 @@ if !isdir(datapath); mkdir(datapath); end nothing # hide # and then create Output. -get_sol(prob) = sol # extracts the Fourier-transformed solution -get_u(prob) = irfft(im * grid.l .* sqrt.(grid.invKrsq) .* sol, grid.nx) +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution +get_u(prob) = irfft(im * prob.grid.l .* sqrt.(prob.grid.invKrsq) .* prob.sol, prob.grid.nx) + out = Output(prob, dataname, (:sol, get_sol), (:u, get_u)) nothing # hide @@ -133,56 +139,42 @@ nothing # hide # We define a function that plots the buoyancy field and the time evolution of kinetic energy # and buoyancy variance. -function plot_output(prob) - b = prob.vars.b - - pb = heatmap(x, y, Array(b'), - aspectratio = 1, - c = :deep, - clim = (0, 1), - xlims = (-grid.Lx/2, grid.Lx/2), - ylims = (-grid.Ly/2, grid.Ly/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "buoyancy bₛ", - framestyle = :box) - - pKE = plot(1, - label = "kinetic energy ∫½(uₛ²+vₛ²)dxdy/L²", - linewidth = 2, - legend = :bottomright, - alpha = 0.7, - xlims = (0, tf), - ylims = (0, 1e-2), - xlabel = "t") - - pb² = plot(1, - label = "buoyancy variance ∫bₛ²dxdy/L²", - linecolor = :red, - legend = :bottomright, - linewidth = 2, - alpha = 0.7, - xlims = (0, tf), - ylims = (0, 2e-2), - xlabel = "t") - - layout = @layout [a{0.5w} Plots.grid(2, 1)] - p = plot(pb, pKE, pb², layout=layout, size = (900, 500)) - - return p -end -nothing # hide +b = Observable(Array(vars.b)) + +ke = Observable([Point2f(KE.t[1], KE.data[1])]) +b² = Observable([Point2f(B.t[1], B.data[1])]) + +title_b = Observable("buoyancy, t=" * @sprintf("%.2f", clock.t)) + +fig = Figure(resolution = (900, 600)) + +axb = Axis(fig[1:2, 1]; + xlabel = "x", + ylabel = "y", + title = title_b, + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axE = Axis(fig[1, 2]; + xlabel = "t", + limits = ((0, tf), (0, 2e-2))) + +heatmap!(axb, x, y, b; + colormap = :deep, colorrange = (0, 1)) + +hE = lines!(axE, ke; linewidth = 3) +hb² = lines!(axE, b²; linewidth = 3) + +Legend(fig[2, 2], [hE, hb²], ["kinetic energy ∫½(uₛ²+vₛ²)dxdy/L²", "buoyancy variance ∫bₛ²dxdy/L²"]) + +fig # ## Time-stepping the `Problem` forward and create animation by updating the plot. startwalltime = time() -p = plot_output(prob) - -anim = @animate for j = 0:round(Int, nsteps/nsubs) +record(fig, "sqg_ellipticalvortex.mp4", 0:round(Int, nsteps/nsubs), framerate = 14) do j if j % (500 / nsubs) == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) @@ -197,65 +189,55 @@ anim = @animate for j = 0:round(Int, nsteps/nsubs) println(log2) end - p[1][1][:z] = Array(vars.b) - p[1][:title] = "buoyancy, t=" * @sprintf("%.2f", clock.t) - push!(p[2][1], KE.t[KE.i], KE.data[KE.i]) - push!(p[3][1], B.t[B.i], B.data[B.i]) + b[] = vars.b + + ke[] = push!(ke[], Point2f(KE.t[KE.i], KE.data[KE.i])) + b²[] = push!(b²[], Point2f(B.t[B.i], B.data[B.i])) + + title_b[] = "buoyancy, t=" * @sprintf("%.2f", clock.t) stepforward!(prob, diags, nsubs) SurfaceQG.updatevars!(prob) end +nothing # hide + +# ![](sqg_ellipticalvortex.mp4) -mp4(anim, "sqg_ellipticalvortex.mp4", fps=14) # Let's see how all flow fields look like at the end of the simulation. -pu = heatmap(x, y, Array(vars.u'), - aspectratio = 1, - c = :balance, - clim = (-maximum(abs.(vars.u)), maximum(abs.(vars.u))), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", - framestyle = :box) - -pv = heatmap(x, y, Array(vars.v'), - aspectratio = 1, - c = :balance, - clim = (-maximum(abs.(vars.v)), maximum(abs.(vars.v))), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", - framestyle = :box) - -pb = heatmap(x, y, Array(vars.b'), - aspectratio = 1, - c = :deep, - clim = (0, 1), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "bₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", - framestyle = :box) +fig = Figure(resolution = (800, 380)) + +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axb = Axis(fig[1, 1]; title = "bₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", axis_kwargs...) +axu = Axis(fig[1, 2]; title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", axis_kwargs...) +axv = Axis(fig[1, 3]; title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", axis_kwargs...) + +hb = heatmap!(axb, x, y, Array(vars.b); + colormap = :deep, colorrange = (0, 1)) + +Colorbar(fig[2, 1], hb, vertical = false) + +hu = heatmap!(axu, x, y, Array(vars.u); + colormap = :balance, colorrange = (-maximum(abs.(vars.u)), maximum(abs.(vars.u)))) + +Colorbar(fig[2, 2], hu, vertical = false) + +hv = heatmap!(axv, x, y, Array(vars.v); + colormap = :balance, colorrange = (-maximum(abs.(vars.v)), maximum(abs.(vars.v)))) + +Colorbar(fig[2, 3], hv, vertical = false) -layout = @layout [a{0.5h}; b{0.5w} c{0.5w}] +fig -plot_final = plot(pb, pu, pv, layout=layout, size = (800, 800)) # ## Save -# Last we can save the output by calling +# We can save the output at the end of the simulation by calling # ```julia -# saveoutput(out)` +# saveoutput(out) # ``` diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index b5c66b74..2ad6efea 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -1,4 +1,4 @@ -# # 2D decaying turbulence +# # [2D decaying turbulence](@id twodnavierstokes_decaying_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/twodnavierstokes_decaying.ipynb). # @@ -10,13 +10,13 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Printf, Random, Plots" +# pkg"add GeophysicalFlows, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, Printf, Random, Plots +using GeophysicalFlows, Printf, Random, CairoMakie using Random: seed! using GeophysicalFlows: peakedisotropicspectrum @@ -48,12 +48,13 @@ nothing # hide # at every time-step that removes enstrophy at high wavenumbers and, thereby, # stabilize the problem, despite that we use the default viscosity coefficient `ν=0`. -prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ny=n, Ly=L, dt=dt, stepper="FilteredRK4") +prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ny=n, Ly=L, dt, stepper="FilteredRK4") nothing # hide # Next we define some shortcuts for convenience. sol, clock, vars, grid = prob.sol, prob.clock, prob.vars, prob.grid -x, y = grid.x, grid.y +x, y = grid.x, grid.y +Lx, Ly = grid.Lx, grid.Ly nothing # hide @@ -69,25 +70,26 @@ nothing # hide # Let's plot the initial vorticity field. Note that when plotting, we decorate the variable # to be plotted with `Array()` to make sure it is brought back on the CPU when `vars` live on # the GPU. -heatmap(x, y, Array(vars.ζ'), - 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", + +fig = Figure() +ax = Axis(fig[1, 1]; + xlabel = "x", + ylabel = "y", title = "initial vorticity", - framestyle = :box) + aspect = 1, + limits = ((-L/2, L/2), (-L/2, L/2))) + +heatmap!(ax, x, y, Array(vars.ζ'); + colormap = :balance, colorrange = (-40, 40)) + +fig # ## Diagnostics # Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top. -E = Diagnostic(TwoDNavierStokes.energy, prob; nsteps=nsteps) -Z = Diagnostic(TwoDNavierStokes.enstrophy, prob; nsteps=nsteps) +E = Diagnostic(TwoDNavierStokes.energy, prob; nsteps) +Z = Diagnostic(TwoDNavierStokes.enstrophy, prob; nsteps) diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -108,7 +110,8 @@ nothing # hide # And then create Output get_sol(prob) = prob.sol # extracts the Fourier-transformed solution -get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx) +get_u(prob) = irfft(im * prob.grid.l .* prob.grid.invKrsq .* prob.sol, prob.grid.nx) + out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) saveproblem(out) nothing # hide @@ -119,31 +122,33 @@ nothing # hide # We initialize a plot with the vorticity field and the time-series of # energy and enstrophy diagnostics. -p1 = heatmap(x, y, Array(vars.ζ'), - 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)/E(0)" "enstrophy Z(t)/Z(0)"], - legend = :right, - linewidth = 2, - alpha = 0.7, - xlabel = "t", - xlims = (0, 41), - ylims = (0, 1.1)) - -l = @layout Plots.grid(1, 2) -p = plot(p1, p2, layout = l, size = (800, 360)) +ζ = Observable(Array(vars.ζ)) +title_ζ = Observable("vorticity, t=" * @sprintf("%.2f", clock.t)) + +energy = Observable(Point2f[(E.t[1], E.data[1] / E.data[1])]) +enstrophy = Observable(Point2f[(Z.t[1], Z.data[1] / Z.data[1])]) + +fig = Figure(resolution = (800, 360)) + +axζ = Axis(fig[1, 1]; + xlabel = "x", + ylabel = "y", + title = title_ζ, + aspect = 1, + limits = ((-L/2, L/2), (-L/2, L/2))) + +ax2 = Axis(fig[1, 2], + xlabel = "t", + limits = ((-0.5, 40.5), (0, 1.05))) + +heatmap!(axζ, x, y, ζ; + colormap = :balance, colorrange = (-40, 40)) +hE = lines!(ax2, energy; linewidth = 3) +hZ = lines!(ax2, enstrophy; linewidth = 3, color = :red) +Legend(fig[1, 3], [hE, hZ], ["E(t)/E(0)", "Z(t)/Z(0)"]) + +fig # ## Time-stepping the `Problem` forward @@ -151,7 +156,7 @@ p = plot(p1, p2, layout = l, size = (800, 360)) startwalltime = time() -anim = @animate for j = 0:Int(nsteps/nsubs) +record(fig, "twodturb.mp4", 0:Int(nsteps/nsubs), framerate = 18) do j if j % (1000 / nsubs) == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) @@ -161,22 +166,19 @@ anim = @animate for j = 0:Int(nsteps/nsubs) println(log) end - p[1][1][:z] = Array(vars.ζ) - p[1][:title] = "vorticity, t=" * @sprintf("%.2f", clock.t) - push!(p[2][1], E.t[E.i], E.data[E.i]/E.data[1]) - push!(p[2][2], Z.t[Z.i], Z.data[Z.i]/Z.data[1]) + ζ[] = vars.ζ + + energy[] = push!(energy[], Point2f(E.t[E.i], E.data[E.i] / E.data[1])) + enstrophy[] = push!(enstrophy[], Point2f(Z.t[E.i], Z.data[Z.i] / Z.data[1])) + + title_ζ[] = "vorticity, t=" * @sprintf("%.2f", clock.t) stepforward!(prob, diags, nsubs) TwoDNavierStokes.updatevars!(prob) end +nothing # hide -mp4(anim, "twodturb.mp4", fps=18) - - -# Last we can save the output by calling -# ```julia -# saveoutput(out)` -# ``` +# ![](twodturb.mp4) # ## Radial energy spectrum @@ -184,17 +186,19 @@ mp4(anim, "twodturb.mp4", fps=18) # After the simulation is done we plot the instantaneous radial energy spectrum to illustrate # how `FourierFlows.radialspectrum` can be used, -E = @. 0.5 * (vars.u^2 + vars.v^2) # energy density -Eh = rfft(E) # Fourier transform of energy density -kr, Ehr = FourierFlows.radialspectrum(Eh, grid, refinement=1) # compute radial specturm of `Eh` +E = @. 0.5 * (vars.u^2 + vars.v^2) # energy density +Eh = rfft(E) # Fourier transform of energy density + +## compute radial specturm of `Eh` +kr, Ehr = FourierFlows.radialspectrum(Eh, grid, refinement = 1) nothing # hide # and we plot it. -plot(kr, abs.(Ehr), - linewidth = 2, - alpha = 0.7, - xlabel = "kᵣ", ylabel = "∫ |Ê| kᵣ dk_θ", - xlims = (5e-1, grid.nx), - xscale = :log10, yscale = :log10, - title = "Radial energy spectrum", - legend = false) +lines(kr, vec(abs.(Ehr)); + linewidth = 2, + axis = (xlabel = L"k_r", + ylabel = L"\int |\hat{E}| k_r \mathrm{d}k_\theta", + xscale = log10, + yscale = log10, + title = "Radial energy spectrum", + limits = ((0.3, 1e2), (1e0, 1e5)))) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 73788a4d..b0fc5abd 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -1,4 +1,4 @@ -# # 2D forced-dissipative turbulence +# # [2D forced-dissipative turbulence](@id twodnavierstokes_stochasticforcing_example) # #md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/twodnavierstokes_stochasticforcing.ipynb). # @@ -12,14 +12,16 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, CUDA, Random, Printf, Plots" +# pkg"add GeophysicalFlows, CUDA, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, CUDA, Random, Printf, Plots +using GeophysicalFlows, CUDA, Random, Printf, CairoMakie + parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU @@ -86,7 +88,7 @@ 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", +prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν, nν, μ, nμ, dt, stepper="ETDRK4", calcF=calcF!, stochastic=true) nothing # hide @@ -105,18 +107,19 @@ nothing # hide # it is brought back on the CPU when the variable lives on the GPU. calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) -heatmap(x, y, Array(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, +fig = Figure() + +ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", - title = "a forcing realization", - framestyle = :box) + aspect = 1, + title = "a forcing realization", + limits = ((-L/2, L/2), (-L/2, L/2))) + +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); + colormap = :balance, colorrange = (-200, 200)) + +fig # ## Setting initial conditions @@ -128,8 +131,8 @@ TwoDNavierStokes.set_ζ!(prob, ArrayType(dev)(zeros(grid.nx, grid.ny))) # ## Diagnostics # Create Diagnostics; the diagnostics are aimed to probe the energy budget. -E = Diagnostic(TwoDNavierStokes.energy, prob, nsteps=nsteps) # energy -Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nsteps) # enstrophy +E = Diagnostic(TwoDNavierStokes.energy, prob; nsteps) # energy +Z = Diagnostic(TwoDNavierStokes.enstrophy, prob; nsteps) # enstrophy diags = [E, Z] # a list of Diagnostics passed to `stepforward!` will be updated every timestep. nothing # hide @@ -137,42 +140,45 @@ 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, Array(vars.ζ'), - 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.55)) - -l = @layout Plots.grid(1, 2) -p = plot(p1, p2, layout = l, size = (900, 420)) +# energy and enstrophy diagnostics. To plot energy and enstrophy on the +# same axes we scale enstrophy with ``k_f^2``. + +ζ = Observable(Array(vars.ζ)) +title_ζ = Observable("vorticity, μ t=" * @sprintf("%.2f", μ * clock.t)) + +energy = Observable(Point2f[(μ * E.t[1], E.data[1])]) +enstrophy = Observable(Point2f[(μ * Z.t[1], Z.data[1] / forcing_wavenumber^2)]) + +fig = Figure(resolution = (800, 360)) + +axζ = Axis(fig[1, 1]; + xlabel = "x", + ylabel = "y", + title = title_ζ, + aspect = 1, + limits = ((-L/2, L/2), (-L/2, L/2))) + +ax2 = Axis(fig[1, 2], + xlabel = "μ t", + limits = ((0, 1.1 * μ * nsteps * dt), (0, 0.55))) + +heatmap!(axζ, x, y, ζ; + colormap = :balance, colorrange = (-40, 40)) + +hE = lines!(ax2, energy; linewidth = 3) +hZ = lines!(ax2, enstrophy; linewidth = 3, color = :red) +Legend(fig[1, 3], [hE, hZ], ["energy E(t)" "enstrophy Z(t) / k_f²"]) + +fig # ## Time-stepping the `Problem` forward -# Finally, we time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() -anim = @animate for j = 0:round(Int, nsteps / nsubs) +record(fig, "twodturb_forced.mp4", 0:round(Int, nsteps / nsubs), framerate = 18) do j if j % (1000/nsubs) == 0 cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy]) @@ -181,13 +187,16 @@ anim = @animate for j = 0:round(Int, nsteps / nsubs) println(log) end - p[1][1][:z] = Array(vars.ζ) - 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) + ζ[] = vars.ζ + + energy[] = push!(energy[], Point2f(μ * E.t[E.i], E.data[E.i])) + enstrophy[] = push!(enstrophy[], Point2f(μ * Z.t[E.i], Z.data[Z.i] / forcing_wavenumber^2)) + + title_ζ[] = "vorticity, μ t=" * @sprintf("%.2f", μ * clock.t) stepforward!(prob, diags, nsubs) TwoDNavierStokes.updatevars!(prob) end +nothing # hide -mp4(anim, "twodturb_forced.mp4", fps=18) +# ![](twodturb_forced.mp4) diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 1bdccf72..7707f475 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -1,4 +1,4 @@ -# # 2D forced-dissipative turbulence budgets +# # [2D forced-dissipative turbulence budgets](@id twodnavierstokes_stochasticforcing_budgets_example) # #md # This example can viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/twodnavierstokes_stochasticforcing_budgets.ipynb). # @@ -14,14 +14,16 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, CUDA, Random, Printf, Plots" +# pkg"add GeophysicalFlows, CUDA, Random, Printf, CairoMakie" # ``` # ## Let's begin -# Let's load `GeophysicalFlows.jl` and some other needed packages. +# Let's load `GeophysicalFlows.jl` and some other packages we need. # -using GeophysicalFlows, CUDA, Random, Printf, Plots +using GeophysicalFlows, CUDA, Random, Printf, CairoMakie + parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU @@ -88,14 +90,15 @@ 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", +prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν, nν, μ, nμ, 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 -x, y = grid.x, grid.y +x, y = grid.x, grid.y +Lx, Ly = grid.Lx, grid.Ly nothing # hide @@ -107,18 +110,19 @@ nothing # hide # it is brought back on the CPU when the variable lives on the GPU. calcF!(vars.Fh, sol, 0.0, clock, vars, params, grid) -heatmap(x, y, Array(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, +fig = Figure() + +ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", - title = "a forcing realization", - framestyle = :box) + aspect = 1, + title = "a forcing realization", + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); + colormap = :balance, colorrange = (-200, 200)) + +fig # ## Setting initial conditions @@ -141,114 +145,9 @@ Wᶻ = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, n 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 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) - 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) - - pζ = heatmap(x, y, Array(vars.ζ'), - 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(pζ, 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. +# We step the `Problem` forward in time. startwalltime = time() for i = 1:ns @@ -266,13 +165,102 @@ end # ## Plot -# And now let's see what we got. First we plot the final snapshot -# of the vorticity field. +# Now let's see the final snapshot of the vorticity. + +fig = Figure(resolution = (400, 400)) + +ax = Axis(fig[1, 1]; + xlabel = "x", + ylabel = "y", + title = "∇²ψ(x, y, μt=" * @sprintf("%.2f", μ * clock.t) * ")", + aspect = 1, + limits = ((-L/2, L/2), (-L/2, L/2))) + +heatmap!(ax, x, y, Array(vars.ζ); + colormap = :viridis, colorrange = (-25, 25)) + +fig + +# And finally, we plot the evolution of the energy and enstrophy diagnostics 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``. + +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) + +t = E.t[2:E.i] + +fig = Figure(resolution = (800, 1100)) + +axis_kwargs = (xlabel = "μ t", ) + +ax1E = Axis(fig[1, 1]; ylabel = "energy sources/sinks", axis_kwargs...) +ax2E = Axis(fig[3, 1]; ylabel = "dE/dt", axis_kwargs...) +ax3E = Axis(fig[5, 1]; axis_kwargs...) + +ax1Z = Axis(fig[1, 2]; axis_kwargs...) +ax2Z = Axis(fig[3, 2]; axis_kwargs...) +ax3Z = Axis(fig[5, 2]; axis_kwargs...) + +hWᵋ = lines!(ax1E, t, Wᵋ[2:E.i]; linestyle = :solid) +hε = lines!(ax1E, t, ε .+ 0t; linestyle = :dash) +hDᵋ = lines!(ax1E, t, Dᵋ[1:E.i-1]; linestyle = :solid) +hRᵋ = lines!(ax1E, t, Rᵋ[1:E.i-1]; linestyle = :solid) + +Legend(fig[2, 1], + [hWᵋ, hε, hDᵋ, hRᵋ], + ["energy work, Wᵋ" "ensemble mean energy work, " "dissipation, Dᵋ" "drag, Rᵋ = - 2μE"]) + +hc = lines!(ax2E, t, dEdt_computed; linestyle = :solid) +hn = lines!(ax2E, t, dEdt_numerical; linestyle = :dash) + +Legend(fig[4, 1], + [hc, hn], + ["computed Wᵋ-Dᵋ" "numerical dE/dt"]) + +hr = lines!(ax3E, t, residual_E) + +Legend(fig[6, 1], + [hr], + ["residual"]) + +hWᶻ = lines!(ax1Z, t, Wᶻ[2:Z.i]; linestyle = :solid) +hεᶻ = lines!(ax1Z, t, εᶻ .+ 0t; linestyle = :dash) +hDᶻ = lines!(ax1Z, t, Dᶻ[1:Z.i-1]; linestyle = :solid) +hRᶻ = lines!(ax1Z, t, Rᶻ[1:Z.i-1]; linestyle = :solid) + +Legend(fig[2, 2], + [hWᶻ, hεᶻ, hDᶻ, hRᶻ], + ["enstrophy work, Wᶻ" "ensemble mean enstophy work, " "dissipation, Dᶻ" "drag, Rᶻ = - 2μZ"]) + +hcᶻ = lines!(ax2Z, t, dZdt_computed; linestyle = :solid) +hnᶻ = lines!(ax2Z, t, dZdt_numerical; linestyle = :dash) -pζ, pbudgets = computetendencies_and_makeplot(prob, diags) +Legend(fig[4, 2], + [hcᶻ, hnᶻ], + ["computed Wᶻ-Dᶻ" "numerical dZ/dt"]) -pζ +hrᶻ = lines!(ax3Z, t, residual_Z) -# And finaly the energy and enstrophy budgets. +Legend(fig[6, 2], + [hr], + ["residual"]) -pbudgets +fig