From 90be2beb1e67121c0b0a0dd1aa92f878dd10e504 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 26 May 2022 09:16:04 +1000 Subject: [PATCH 01/44] use FourierFlows v0.9.1 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f46016a8..c9a15455 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" FFTW = "^1" -FourierFlows = "^0.9" +FourierFlows = "^0.9.1" JLD2 = "^0.1, ^0.2, ^0.3, ^0.4" Reexport = "^0.2, ^1" SpecialFunctions = "^0.10, ^1, 2" From 1d6943b2a44697ebcddca150113bbbca29a7f175 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 26 May 2022 09:16:25 +1000 Subject: [PATCH 02/44] minor fig size adjustment --- examples/twodnavierstokes_stochasticforcing_budgets.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 1bdccf72..69c10298 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -237,7 +237,7 @@ function computetendencies_and_makeplot(prob, diags) layout = @layout Plots.grid(3, 2) - pbudgets = plot(p1E, p1Z, p2E, p2Z, p3E, p3Z, layout=layout, size = (900, 1200)) + pbudgets = plot(p1E, p1Z, p2E, p2Z, p3E, p3Z, layout=layout, size = (800, 1100)) return pζ, pbudgets end @@ -273,6 +273,6 @@ pζ, pbudgets = computetendencies_and_makeplot(prob, diags) pζ -# And finaly the energy and enstrophy budgets. +# And finaly, we plot the energy and enstrophy budgets. pbudgets From 8ff164022636215becc96bd397226f4b919ef0b3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 12 Jul 2022 18:12:14 +0300 Subject: [PATCH 03/44] use FourierFlows v0.9.2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d8c0f1ee..0e68b3ab 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.1" +FourierFlows = "^0.9.2" JLD2 = "^0.1, ^0.2, ^0.3, ^0.4" Reexport = "^0.2, ^1" SpecialFunctions = "^0.10, ^1, 2" From 4f0d00d51f8c7153e7a8cf1e61010f0b52e4c7b6 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 27 Aug 2022 16:35:11 +1000 Subject: [PATCH 04/44] use Makie for plots/animations --- examples/barotropicqgql_betaforced.jl | 236 +++++++++++++------------- 1 file changed, 118 insertions(+), 118 deletions(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 8259fbba..730007cc 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -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,13 +98,14 @@ 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 @@ -110,18 +114,19 @@ 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, irfft(vars.Fh, grid.nx); + colormap = :balance, colorrange = (-8, 8)) + +fig # ## Setting initial conditions @@ -133,17 +138,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 +168,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 +188,87 @@ 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 -nothing # hide +title_ζ = Observable(@sprintf("vorticity, μt = %.2f", μ * clock.t)) +title_ψ = "streamfunction ψ" + +fig = Figure(resolution=(1000, 600)) + +axζ = Axis(fig[1, 1], + xlabel = "x", + ylabel = "y", + aspect = 1, + title = title_ζ, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axψ = Axis(fig[2, 1], + xlabel = "x", + ylabel = "y", + aspect = 1, + title = title_ψ, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +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(@. ζ̄ + ζ′) +ψ̄, ψ′= prob.vars.Psi, prob.vars.psi +ψ = Observable(@. ψ̄ + ψ′) +ζ̄ₘ = Observable(vec(mean(ζ̄, dims=1))) +ūₘ = Observable(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) + +fig # ## Time-stepping the `Problem` forward # 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, "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,20 +279,21 @@ 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 -mp4(anim, "barotropicqgql_betaforced.mp4", fps=18) - # ## Save From 2da74775e28f0f14edb320f5df3f6471ed9b32a1 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Aug 2022 08:03:46 +1000 Subject: [PATCH 05/44] cosmetics --- examples/multilayerqg_2layer.jl | 25 ++++++++++--------- examples/singlelayerqg_betadecay.jl | 3 +-- examples/singlelayerqg_betaforced.jl | 2 +- ...ecaying_barotropic_equivalentbarotropic.jl | 6 ++--- examples/singlelayerqg_decaying_topography.jl | 2 +- examples/surfaceqg_decaying.jl | 2 +- examples/twodnavierstokes_decaying.jl | 2 +- .../twodnavierstokes_stochasticforcing.jl | 6 ++--- ...dnavierstokes_stochasticforcing_budgets.jl | 2 +- 9 files changed, 25 insertions(+), 25 deletions(-) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 644c8921..f066966b 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -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. # -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 diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index ddbfd7e2..790c0a86 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -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 diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 7e19e26e..3928a574 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -94,7 +94,7 @@ 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 diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index b2e01ca2..febf661f 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -55,9 +55,9 @@ nothing # hide # 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) + dt, stepper="FilteredRK4", aliased_fraction=0) +prob_eqbqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, deformation_radius, + dt, stepper="FilteredRK4", aliased_fraction=0) nothing # hide diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 1de9b047..25e95437 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -57,7 +57,7 @@ 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 diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index 492f140a..fb655d24 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -58,7 +58,7 @@ 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. diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index b5c66b74..236062d3 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -48,7 +48,7 @@ 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. diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 73788a4d..6dfeb148 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -86,7 +86,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 @@ -128,8 +128,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 diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 69c10298..6ef28221 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -88,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 From 63e893b6c892fa6130c31d315d828469c1824742 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Aug 2022 08:06:18 +1000 Subject: [PATCH 06/44] add CairoMakie --- docs/Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index 0511b5cd..4dd4f832 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,8 @@ [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" From 0b2b5e3837b47119bc3fdc34aa486aacfa4a4600 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Aug 2022 15:50:31 +1000 Subject: [PATCH 07/44] bump patch release --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 0e68b3ab..d8e56b02 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "GeophysicalFlows" uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e" license = "MIT" authors = ["Navid C. Constantinou ", "Gregory L. Wagner ", "and co-contributors"] -version = "0.14.1" +version = "0.14.2" [deps] CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" From 23b56c3688eb42ee238380c00e8409efdbfaa3a4 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Aug 2022 15:50:53 +1000 Subject: [PATCH 08/44] use Makie for plots/animations --- examples/multilayerqg_2layer.jl | 177 +++++++++++++++++++------------- 1 file changed, 103 insertions(+), 74 deletions(-) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index f066966b..6825fb7c 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -135,104 +135,133 @@ nothing # hide # 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(vars.q[:, :, 1]) +ψ₁ = Observable(vars.ψ[:, :, 1]) +q₂ = Observable(vars.q[:, :, 2]) +ψ₂ = Observable(vars.ψ[:, :, 2]) +function compute_levels(maxf) + levelsf = @lift collect(range(-$maxf, stop = $maxf, length=8)) + levelsf⁺ = @lift collect(range($maxf/7, stop = $maxf, length=4)) + levelsf⁻ = @lift collect(range(-$maxf, stop = -$maxf/7, length=4)) + + 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)) + +axq₁ = Axis(fig[1, 1], + xlabel = "x", + ylabel = "y", + title = "q₁", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axψ₁ = Axis(fig[2, 1], + xlabel = "x", + ylabel = "y", + title = "ψ₁", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axq₂ = Axis(fig[1, 2], + xlabel = "x", + ylabel = "y", + title = "q₂", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axψ₂ = Axis(fig[2, 2], + xlabel = "x", + ylabel = "y", + aspect = 1, + title = "ψ₂", + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +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 -mp4(anim, "multilayerqg_2layer.mp4", fps=18) - # ## Save # Finally, we can save, e.g., the last snapshot via From 3a46a5c29984c9a5ef6578763e34aad21a3be691 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Aug 2022 16:12:27 +1000 Subject: [PATCH 09/44] use Makie for plots/animations --- examples/singlelayerqg_betadecay.jl | 200 +++++++++++++--------------- 1 file changed, 93 insertions(+), 107 deletions(-) diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 790c0a86..7541be69 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -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. # -using GeophysicalFlows, Plots, Printf, Random +using GeophysicalFlows, CairoMakie, Printf, Random using Statistics: mean @@ -86,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 @@ -141,7 +134,7 @@ 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 @@ -151,67 +144,63 @@ nothing # hide # 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 +Lx, Ly = grid.Lx, grid.Ly + +title_q = Observable(@sprintf("vorticity, t = %.2f", clock.t)) +title_ψ = "streamfunction ψ" + +fig = Figure(resolution=(800, 720)) + +axq = Axis(fig[1, 1], + xlabel = "x", + ylabel = "y", + aspect = 1, + title = title_q, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +axψ = Axis(fig[2, 1], + xlabel = "x", + ylabel = "y", + aspect = 1, + title = title_ψ, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) + +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(vars.q) +ψ = Observable(vars.ψ) +q̄ₘ = Observable(vec(mean(vars.q, dims=1))) +ūₘ = Observable(vec(mean(vars.u, dims=1))) + +heatmap!(axq, x, y, q; + colormap = :balance, colorrange = (-12, 12)) + +contourf!(axψ, x, y, ψ; + levels = collect(range(-0.7, stop=0.7, length=20)), + colormap = :viridis, + colorrange = (-0.35, 0.35)) + +contour!(axψ, x, y, ψ; + levels = collect(range(-0.7, stop=0.7, length=20)), + 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 @@ -220,10 +209,9 @@ nothing # hide 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]) @@ -231,21 +219,19 @@ 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 -mp4(anim, "singlelayerqg_betadecay.mp4", fps=8) - # ## Save # Finally, we can save, e.g., the last snapshot via From 6d7c2916cfd8cafb0f30a571f1b1b0342d826b46 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 28 Aug 2022 16:36:28 +1000 Subject: [PATCH 10/44] minor cosmetics --- examples/singlelayerqg_betaforced.jl | 8 ++++---- examples/singlelayerqg_decaying_topography.jl | 6 +++--- examples/surfaceqg_decaying.jl | 4 ++-- examples/twodnavierstokes_decaying.jl | 4 ++-- examples/twodnavierstokes_stochasticforcing.jl | 4 ++-- 5 files changed, 13 insertions(+), 13 deletions(-) diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 3928a574..83c46efc 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -132,8 +132,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) +Z = Diagnostic(SingleLayerQG.enstrophy, prob; nsteps) diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -153,8 +153,8 @@ 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) +get_sol(prob) = prob.sol # extracts the Fourier-transformed solution +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)) nothing # hide diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 25e95437..6a3bdaa7 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -140,8 +140,8 @@ p = plot(p1, p2, layout=layout, size = (800, 360)) # ## 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,7 +158,7 @@ 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 diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index fb655d24..c6a32cf4 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -122,8 +122,8 @@ 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 diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index 236062d3..c0345d1c 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -86,8 +86,8 @@ heatmap(x, y, Array(vars.ζ'), # ## 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 diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 6dfeb148..237e393f 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -128,8 +128,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) # energy -Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, 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 From d8e09700b1160f7d803f68cd9304c7c251267c2f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 07:26:37 +1000 Subject: [PATCH 11/44] show mp4 --- examples/barotropicqgql_betaforced.jl | 4 +++- examples/multilayerqg_2layer.jl | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 730007cc..24c6d097 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -257,7 +257,7 @@ lines!(axū, 0y, y; linewidth = 1, linestyle=:dash) lines!(axE, μt, energy; linewidth = 3) lines!(axZ, μt, enstrophy; linewidth = 3, color = :red) -fig +nothing # hide # ## Time-stepping the `Problem` forward @@ -294,6 +294,8 @@ record(fig, "barotropicqgql_betaforced.mp4", frames, framerate = 18) do j BarotropicQGQL.updatevars!(prob) end +# ![](barotropicqgql_betaforced.mp4) + # ## Save diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 6825fb7c..ec873c1d 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -262,6 +262,8 @@ record(fig, "multilayerqg_2layer.mp4", frames, framerate = 18) do j MultiLayerQG.updatevars!(prob) end +# ![](multilayerqg_2layer.mp4) + # ## Save # Finally, we can save, e.g., the last snapshot via From b44e172ed5bfa2dae065694d5f5fd929b00b07b9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 07:44:36 +1000 Subject: [PATCH 12/44] show mp4 --- examples/singlelayerqg_betadecay.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 7541be69..8c2e3512 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -232,6 +232,9 @@ record(fig, "singlelayerqg_betadecay.mp4", frames, framerate = 8) do j SingleLayerQG.updatevars!(prob) end +# ![](singlelayerqg_betadecay.mp4) + + # ## Save # Finally, we can save, e.g., the last snapshot via From c76700dcd971af19727339aa41cc0915a1143848 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 07:44:56 +1000 Subject: [PATCH 13/44] cleanup --- examples/multilayerqg_2layer.jl | 38 ++++++++++++++++++++++----------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index ec873c1d..7fe79222 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -131,8 +131,9 @@ 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. Lx, Ly = grid.Lx, grid.Ly @@ -144,10 +145,15 @@ q₁ = Observable(vars.q[:, :, 1]) q₂ = Observable(vars.q[:, :, 2]) ψ₂ = Observable(vars.ψ[:, :, 2]) -function compute_levels(maxf) - levelsf = @lift collect(range(-$maxf, stop = $maxf, length=8)) - levelsf⁺ = @lift collect(range($maxf/7, stop = $maxf, length=4)) - levelsf⁻ = @lift collect(range(-$maxf, stop = -$maxf/7, length=4)) +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=nlevels/2)) + + # only negative + levelsf⁻ = @lift collect(range(-$maxf, stop = -$maxf/(nlevels-1), length=nlevels/2)) return levelsf, levelsf⁺, levelsf⁻ end @@ -209,13 +215,19 @@ 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) -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) @@ -256,7 +268,7 @@ record(fig, "multilayerqg_2layer.mp4", frames, framerate = 18) do j 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) + title_KE[] = @sprintf("μ t = %.2f", μ * clock.t) stepforward!(prob, diags, nsubs) MultiLayerQG.updatevars!(prob) From 2ffd6c4bc377f22e5922af2c643508027aa6085a Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 08:24:38 +1000 Subject: [PATCH 14/44] cleanup axis kwargs --- examples/multilayerqg_2layer.jl | 39 ++++++++++----------------------- 1 file changed, 12 insertions(+), 27 deletions(-) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 7fe79222..555fd241 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -170,33 +170,18 @@ PE = Observable(Point2f[(μ * E.t[1], E.data[1][2])]) fig = Figure(resolution=(1000, 600)) -axq₁ = Axis(fig[1, 1], - xlabel = "x", - ylabel = "y", - title = "q₁", - aspect = 1, - limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) - -axψ₁ = Axis(fig[2, 1], - xlabel = "x", - ylabel = "y", - title = "ψ₁", - aspect = 1, - limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) - -axq₂ = Axis(fig[1, 2], - xlabel = "x", - ylabel = "y", - title = "q₂", - aspect = 1, - limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) - -axψ₂ = Axis(fig[2, 2], - xlabel = "x", - ylabel = "y", - aspect = 1, - title = "ψ₂", - limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) +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", From 6d75bda26be113dc7f05831a43c42a5d9409cafe Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 08:26:21 +1000 Subject: [PATCH 15/44] remove package from docs/Project.toml --- docs/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 4dd4f832..edddbe2b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,7 +2,6 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" From fdad6bad94fa709e99e5fa1a0cd5dc17b8a9f07e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 09:58:18 +1000 Subject: [PATCH 16/44] some fixes --- examples/multilayerqg_2layer.jl | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 555fd241..566a60d2 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -140,20 +140,20 @@ Lx, Ly = grid.Lx, grid.Ly title_KE = Observable(@sprintf("μt = %.2f", μ * clock.t)) -q₁ = Observable(vars.q[:, :, 1]) -ψ₁ = Observable(vars.ψ[:, :, 1]) -q₂ = Observable(vars.q[:, :, 2]) -ψ₂ = Observable(vars.ψ[:, :, 2]) +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=nlevels/2)) + 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=nlevels/2)) + levelsf⁻ = @lift collect(range(-$maxf, stop = -$maxf/(nlevels-1), length=Int(nlevels/2))) return levelsf, levelsf⁺, levelsf⁻ end @@ -175,13 +175,13 @@ axis_kwargs = (xlabel = "x", aspect = 1, limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) -axq₁ = Axis(fig[1, 1], title = "q₁", axis_kwargs...) +axq₁ = Axis(fig[1, 1]; title = "q₁", axis_kwargs...) -axψ₁ = Axis(fig[2, 1], title = "ψ₁", axis_kwargs...) +axψ₁ = Axis(fig[2, 1]; title = "ψ₁", axis_kwargs...) -axq₂ = Axis(fig[1, 2], title = "q₂", axis_kwargs...) +axq₂ = Axis(fig[1, 2]; title = "q₂", axis_kwargs...) -axψ₂ = Axis(fig[2, 2], title = "ψ₂", axis_kwargs...) +axψ₂ = Axis(fig[2, 2]; title = "ψ₂", axis_kwargs...) axKE = Axis(fig[1, 3], xlabel = "μ t", @@ -241,10 +241,10 @@ record(fig, "multilayerqg_2layer.mp4", frames, framerate = 18) do j println(log) end - q₁[] = vars.q[:, :, 1] - ψ₁[] = vars.ψ[:, :, 1] - q₂[] = vars.q[:, :, 2] - ψ₂[] = vars.ψ[:, :, 2] + q₁[] = Array(vars.q[:, :, 1]) + ψ₁[] = Array(vars.ψ[:, :, 1]) + q₂[] = Array(vars.q[:, :, 2]) + ψ₂[] = Array(vars.ψ[:, :, 2]) maxψ₁[] = maximum(abs, vars.ψ[:, :, 1]) maxψ₂[] = maximum(abs, vars.ψ[:, :, 2]) From 97eb67943a81d19b3bcde9b8af1a76a5a7a4bdb1 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 10:04:23 +1000 Subject: [PATCH 17/44] some cleanup --- examples/barotropicqgql_betaforced.jl | 19 +++++--------- examples/singlelayerqg_betadecay.jl | 37 +++++++++++---------------- 2 files changed, 22 insertions(+), 34 deletions(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 24c6d097..17926449 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -193,19 +193,14 @@ title_ψ = "streamfunction ψ" fig = Figure(resolution=(1000, 600)) -axζ = Axis(fig[1, 1], - xlabel = "x", - ylabel = "y", - aspect = 1, - title = title_ζ, - limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) -axψ = Axis(fig[2, 1], - xlabel = "x", - ylabel = "y", - aspect = 1, - title = title_ψ, - 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 ζ", diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 8c2e3512..46526014 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -151,19 +151,14 @@ title_ψ = "streamfunction ψ" fig = Figure(resolution=(800, 720)) -axq = Axis(fig[1, 1], - xlabel = "x", - ylabel = "y", - aspect = 1, - title = title_q, - limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) +axis_kwargs = (xlabel = "x", + ylabel = "y", + aspect = 1, + limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) -axψ = Axis(fig[2, 1], - xlabel = "x", - ylabel = "y", - aspect = 1, - title = title_ψ, - 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", @@ -177,28 +172,26 @@ axū = Axis(fig[2, 2], aspect = 1, limits = ((-0.5, 0.5), (-Ly/2, Ly/2))) -q = Observable(vars.q) -ψ = Observable(vars.ψ) +q = Observable(vars.q) +ψ = Observable(vars.ψ) q̄ₘ = Observable(vec(mean(vars.q, dims=1))) ūₘ = Observable(vec(mean(vars.u, dims=1))) heatmap!(axq, x, y, q; colormap = :balance, colorrange = (-12, 12)) -contourf!(axψ, x, y, ψ; - levels = collect(range(-0.7, stop=0.7, length=20)), - colormap = :viridis, - colorrange = (-0.35, 0.35)) +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 = collect(range(-0.7, stop=0.7, length=20)), - color = :black) + levels, color = :black) lines!(axq̄, q̄ₘ, y; linewidth = 3) -lines!(axq̄, 0y, y; linewidth = 1, linestyle=:dash) +lines!(axq̄, 0y, y; linewidth = 1, linestyle = :dash) lines!(axū, ūₘ, y; linewidth = 3) -lines!(axū, 0y, y; lindewidth = 1, linestyle=:dash) +lines!(axū, 0y, y; lindewidth = 1, linestyle = :dash) fig From 8032bf9b5d5ee6157bdcfa0c0d64a76fab482701 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 11:42:53 +1000 Subject: [PATCH 18/44] cleanup --- ...glelayerqg_decaying_barotropic_equivalentbarotropic.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index febf661f..e30844f0 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -54,10 +54,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, stepper="FilteredRK4", aliased_fraction=0) -prob_eqbqg = SingleLayerQG.Problem(dev; nx=n, Lx=L, deformation_radius, - 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 From 97b43b0df8fccfa868f22a8daabe0679e51baeed Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 11:43:56 +1000 Subject: [PATCH 19/44] use Makie for plot + save output and load --- examples/singlelayerqg_betaforced.jl | 255 +++++++++++++++------------ 1 file changed, 139 insertions(+), 116 deletions(-) diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 83c46efc..198e0248 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -11,15 +11,17 @@ # ```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. # -using GeophysicalFlows, CUDA, Random, Printf, Plots +using GeophysicalFlows, CUDA, JLD2, CairoMakie, Random, Printf using Statistics: mean +using LinearAlgebra: ldiv! + parsevalsum = FourierFlows.parsevalsum # ## Choosing a device: CPU or GPU @@ -33,8 +35,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 @@ -121,7 +124,7 @@ heatmap(x, y, Array(irfft(vars.Fh, grid.nx)'), ylabel = "y", title = "a forcing realization", framestyle = :box) - + # ## Setting initial conditions @@ -132,8 +135,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) -Z = Diagnostic(SingleLayerQG.enstrophy, prob; 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 +147,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, @@ -154,98 +157,27 @@ nothing # hide # and then create Output. get_sol(prob) = prob.sol # extracts the Fourier-transformed solution -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)) -nothing # hide +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 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 +185,129 @@ 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`. We can +# now load the JLD2 output and create a time series for fields that interest us. + +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["diags/energy/t"] +E_data = file["diags/energy/data"] +Z_data = file["diags/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(- 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]) + +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))) + +μt = Observable(μ * E_t[1:1]) +energy = Observable(E_data[1:1]) +enstrophy = Observable(Z_data[1: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, μt, energy; linewidth = 3) +lines!(axZ, μt, 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 + μt.val = μ * E_t[1:i] + energy[] = E_data[1:i] + enstrophy[] = Z_data[1:i] +end + +# ![](singlelayerqg_betaforced.mp4) From 3eb3cb4a256aa1f1b031ee48da490c1a546688cf Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 11:44:31 +1000 Subject: [PATCH 20/44] add JLD2 --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index edddbe2b..20565c7b 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,7 @@ 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" From cf0d602626d172acbf175bb060f750d6c697fce8 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 11:59:20 +1000 Subject: [PATCH 21/44] clarify record --- examples/singlelayerqg_betaforced.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 198e0248..c1fa541c 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -23,6 +23,7 @@ using Statistics: mean using LinearAlgebra: ldiv! parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU @@ -32,11 +33,11 @@ nothing # hide # ## Numerical parameters and time-stepping parameters - n = 128 # 2D resolution: n² grid points + n = 128*4 # 2D resolution: n² grid points stepper = "FilteredRK4" # timestepper - dt = 0.05 # timestep - nsteps = 8000 # total number of timesteps - save_substeps = 10 # number of timesteps after which output is saved + dt = 0.05/4 # timestep + nsteps = 8000*4 # total number of timesteps + save_substeps = 10*4 # number of timesteps after which output is saved nothing # hide From 54565fd76db6fd468d6174c381ad31b082eba110 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 29 Aug 2022 22:08:12 +1000 Subject: [PATCH 22/44] up compat for FourierFlows --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index d8e56b02..e6ccc45f 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.3" JLD2 = "^0.1, ^0.2, ^0.3, ^0.4" Reexport = "^0.2, ^1" SpecialFunctions = "^0.10, ^1, 2" From 36a23c516d73284fd4e1974c14b1a8be012e52a6 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 09:24:10 +1000 Subject: [PATCH 23/44] use Makie for plotting --- ...ecaying_barotropic_equivalentbarotropic.jl | 93 ++++++++++--------- 1 file changed, 48 insertions(+), 45 deletions(-) diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index e30844f0..eca928af 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -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. # -using GeophysicalFlows, Printf, Random, Plots +using GeophysicalFlows, Printf, Random, CairoMakie using GeophysicalFlows: peakedisotropicspectrum +using LinearAlgebra: ldiv! using Random: seed! @@ -91,43 +92,45 @@ 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. +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(relativevorticity(prob_bqg)) +ζ_eqbqg = Observable(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 +140,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 +151,16 @@ 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 -mp4(anim, "singlelayerqg_barotropic_equivalentbarotropic.mp4", fps=18) +# ![](singlelayerqg_barotropic_equivalentbarotropic.mp4) From 821b2282349b87693a9f6abd763496e9a54881e2 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 09:48:45 +1000 Subject: [PATCH 24/44] use Makie for plotting --- examples/singlelayerqg_decaying_topography.jl | 181 +++++++----------- 1 file changed, 73 insertions(+), 108 deletions(-) diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 6a3bdaa7..9bf23f8b 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -10,13 +10,13 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, Plots, Printf, Random, Statistics" +# pkg"add GeophysicalFlows, CairoMakie, Printf, Random, Statistics" # ``` # ## Let's begin # Let's load `GeophysicalFlows.jl` and some other needed packages. # -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,30 @@ nothing # hide # # The topophic PV is prescribed via keyword argument `eta`. prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, eta=topographicPV, - dt, 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") +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, params.eta; + levels = collect(-3:0.4:3), colormap = :balance, colorrange = (-3, 3)) + +fig # ## Setting initial conditions @@ -104,37 +103,40 @@ 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(vars.q) +ψ = Observable(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 @@ -165,56 +167,20 @@ 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. We add the topography contours and +# also the time. + +η = prob.params.eta + +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 @@ -224,10 +190,7 @@ nothing # hide 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,16 @@ 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 -mp4(anim, "singlelayerqg_decaying_topography.mp4", fps=12) +# ![](singlelayerqg_decaying_topography.mp4) From 1b3591593ecefae659b5f9e7e29234d59f70e785 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 10:41:08 +1000 Subject: [PATCH 25/44] use Makie for plotting --- examples/twodnavierstokes_decaying.jl | 118 ++++++++++++++------------ 1 file changed, 65 insertions(+), 53 deletions(-) diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index c0345d1c..008ab195 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -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. # -using GeophysicalFlows, Printf, Random, Plots +using GeophysicalFlows, Printf, Random, CairoMakie using Random: seed! using GeophysicalFlows: peakedisotropicspectrum @@ -53,7 +53,8 @@ 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,18 +70,19 @@ 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 @@ -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,37 @@ 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(vars.ζ) +title_ζ = Observable("vorticity, t=" * @sprintf("%.2f", clock.t)) + +t = Observable(E.t[1:1]) +energy = Observable(E.data[1:1] / E.data[1]) +enstrophy = Observable(Z.data[1:1] / Z.data[1]) + +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 +160,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,16 +170,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]) + ζ[] = Array(vars.ζ) + + t.val = E.t[1:E.i] + 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 -mp4(anim, "twodturb.mp4", fps=18) +# ![](twodturb.mp4) # Last we can save the output by calling @@ -190,11 +202,11 @@ kr, Ehr = FourierFlows.radialspectrum(Eh, grid, refinement=1) # compute radial s 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 = "kᵣ", + ylabel = "∫ |Ê| kᵣ dk_θ", + xscale = log10, + yscale = log10, + title = "Radial energy spectrum", + limits = ((1e-1, 1e2), (1e0, 1e4)))) From 0e3e49cef34191ba57db25bca567d05b6c755c68 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 11:15:09 +1000 Subject: [PATCH 26/44] use Makie for plotting --- examples/twodnavierstokes_decaying.jl | 6 +- .../twodnavierstokes_stochasticforcing.jl | 97 ++++---- ...dnavierstokes_stochasticforcing_budgets.jl | 235 +++++++++--------- 3 files changed, 164 insertions(+), 174 deletions(-) diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index 008ab195..65e0ce35 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -125,10 +125,6 @@ nothing # hide ζ = Observable(vars.ζ) title_ζ = Observable("vorticity, t=" * @sprintf("%.2f", clock.t)) -t = Observable(E.t[1:1]) -energy = Observable(E.data[1:1] / E.data[1]) -enstrophy = Observable(Z.data[1:1] / Z.data[1]) - energy = Observable(Point2f[(E.t[1], E.data[1] / E.data[1])]) enstrophy = Observable(Point2f[(Z.t[1], Z.data[1] / Z.data[1])]) @@ -209,4 +205,4 @@ lines(kr, vec(abs.(Ehr)); xscale = log10, yscale = log10, title = "Radial energy spectrum", - limits = ((1e-1, 1e2), (1e0, 1e4)))) + limits = ((0.3, 1e2), (1e0, 1e5)))) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 237e393f..d2156332 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -12,13 +12,14 @@ # ```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. # -using GeophysicalFlows, CUDA, Random, Printf, Plots +using GeophysicalFlows, CUDA, Random, Printf, CairoMakie + parsevalsum = FourierFlows.parsevalsum # ## Choosing a device: CPU or GPU @@ -105,18 +106,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, irfft(vars.Fh, grid.nx); + colormap = :balance, colorrange = (-200, 200)) + +fig # ## Setting initial conditions @@ -137,33 +139,36 @@ 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(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 @@ -172,7 +177,7 @@ p = plot(p1, p2, layout = l, size = (900, 420)) 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 +186,15 @@ 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) + ζ[] = Array(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 -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 6ef28221..b08f5ad5 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -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. # -using GeophysicalFlows, CUDA, Random, Printf, Plots +using GeophysicalFlows, CUDA, Random, Printf, CairoMakie + parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU @@ -107,18 +109,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, irfft(vars.Fh, grid.nx); + colormap = :balance, colorrange = (-200, 200)) + +fig # ## Setting initial conditions @@ -141,114 +144,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 = (800, 1100)) - - return pζ, pbudgets -end -nothing # hide - - - - # ## Time-stepping the `Problem` forward -# Finally, we time-step the `Problem` forward in time. +# We time-step the `Problem` forward in time. startwalltime = time() for i = 1:ns @@ -266,13 +164,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 finaly, 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, we plot the energy and enstrophy budgets. +Legend(fig[6, 2], + [hr], + ["residual"]) -pbudgets +fig From 6e8b516150f79de8e0ec89cf3c6b2a37ead52e2f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 11:15:22 +1000 Subject: [PATCH 27/44] drop Plots --- docs/Project.toml | 2 -- examples/Project.toml | 5 +---- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 20565c7b..e0a67d1e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,10 +4,8 @@ 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/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" From 7b3d35d6ea3217ad679ea64102de1acf63429abf Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 11:42:39 +1000 Subject: [PATCH 28/44] drop Plots --- docs/make.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) 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 ##### From 740cf019f24577c2757e8b8d8d7e9176fc1acf2b Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 15:30:15 +1000 Subject: [PATCH 29/44] add links to modules --- README.md | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 7b211aac..b4248694 100644 --- a/README.md +++ b/README.md @@ -62,14 +62,14 @@ 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 +* [`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`: a multi-layer quasi-geostrophic model over topography and with the ability +* [`MultiLayerQG`](https://fourierflows.github.io/GeophysicalFlowsDocumentation/stable/modules/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. +* [`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 From 87181e3837fd03266c038510416573212780f72d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 15:31:22 +1000 Subject: [PATCH 30/44] fix + convert surfaceqg example to Makie --- examples/barotropicqgql_betaforced.jl | 2 +- examples/multilayerqg_2layer.jl | 2 +- examples/singlelayerqg_betadecay.jl | 5 +- examples/singlelayerqg_betaforced.jl | 58 +++--- ...ecaying_barotropic_equivalentbarotropic.jl | 2 +- examples/singlelayerqg_decaying_topography.jl | 2 +- examples/surfaceqg_decaying.jl | 195 ++++++++---------- examples/twodnavierstokes_decaying.jl | 3 +- .../twodnavierstokes_stochasticforcing.jl | 2 +- ...dnavierstokes_stochasticforcing_budgets.jl | 2 +- 10 files changed, 127 insertions(+), 146 deletions(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 17926449..65b0e64f 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). # diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 566a60d2..30dfe713 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). # diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 46526014..5742fa3a 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). # @@ -141,8 +141,7 @@ nothing # hide # ## Visualizing the simulation -# We define a function that plots the vorticity and streamfunction and -# their corresponding zonal mean structure. +# We plot the vorticity and streamfunction and their corresponding zonal mean structure. Lx, Ly = grid.Lx, grid.Ly diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index c1fa541c..aa89faff 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). # @@ -33,11 +33,11 @@ nothing # hide # ## Numerical parameters and time-stepping parameters - n = 128*4 # 2D resolution: n² grid points + n = 128 # 2D resolution: n² grid points stepper = "FilteredRK4" # timestepper - dt = 0.05/4 # timestep - nsteps = 8000*4 # total number of timesteps - save_substeps = 10*4 # number of timesteps after which output is saved + dt = 0.05 # timestep + nsteps = 8000 # total number of timesteps + save_substeps = 10 # number of timesteps after which output is saved nothing # hide @@ -104,7 +104,8 @@ 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 @@ -113,18 +114,19 @@ 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 @@ -221,9 +223,10 @@ 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["diags/energy/t"] -E_data = file["diags/energy/data"] -Z_data = file["diags/enstrophy/data"] +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"] @@ -240,6 +243,9 @@ ū = @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", @@ -275,10 +281,6 @@ axZ = Axis(fig[2, 3], aspect = 1, limits = ((-0.1, 4.1), (0, 3.1))) -μt = Observable(μ * E_t[1:1]) -energy = Observable(E_data[1:1]) -enstrophy = Observable(Z_data[1:1]) - heatmap!(axq, x, y, q; colormap = :balance, colorrange = (-8, 8)) @@ -295,8 +297,8 @@ lines!(axq̄, 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) +lines!(axE, energy; linewidth = 3) +lines!(axZ, enstrophy; linewidth = 3, color = :red) fig @@ -306,9 +308,9 @@ fig frames = 1:length(t) record(fig, "singlelayerqg_betaforced.mp4", frames, framerate = 18) do i j[] = i - μt.val = μ * E_t[1:i] - energy[] = E_data[1:i] - enstrophy[] = Z_data[1:i] + + energy[] = push!(energy[], Point2f(μ * E_t[i], E_data[i])) + enstrophy[] = push!(enstrophy[], Point2f(μ * Z_t[i], Z_data[i])) end # ![](singlelayerqg_betaforced.mp4) diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index eca928af..c8bba245 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). # diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 9bf23f8b..241f02cd 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). # diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index c6a32cf4..f463eaba 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. # -using GeophysicalFlows, Plots, Printf, Random +using GeophysicalFlows, CairoMakie, Printf, Random using Statistics: mean using Random: seed! @@ -63,7 +63,8 @@ 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 @@ -124,6 +129,7 @@ nothing # hide # and then create Output. 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(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,61 +189,50 @@ 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[] = Array(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 -mp4(anim, "sqg_ellipticalvortex.mp4", fps=14) +# ![](sqg_ellipticalvortex.mp4) + # 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)) + +kwargs_axis = (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) * ")", kwargs_axis...) +axu = Axis(fig[1, 2]; title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", kwargs_axis...) +axv = Axis(fig[1, 3]; title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", kwargs_axis...) + +hb = heatmap!(axb, x, y, vars.b; + colormap = :deep, colorrange = (0, 1)) + +Colorbar(fig[2, 1], hb, vertical = false) + +hu = heatmap!(axu, x, y, 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, 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 diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index 65e0ce35..b5567be2 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). # @@ -168,7 +168,6 @@ record(fig, "twodturb.mp4", 0:Int(nsteps/nsubs), framerate = 18) do j ζ[] = Array(vars.ζ) - t.val = E.t[1:E.i] 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])) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index d2156332..3c01b64d 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). # diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index b08f5ad5..66920a4b 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). # From 4e7ab38df308f947743a17cb684bb81e98e92791 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 16:09:54 +1000 Subject: [PATCH 31/44] clarify record --- docs/Project.toml | 1 + examples/twodnavierstokes_stochasticforcing.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index e0a67d1e..97319fb9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,7 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 3c01b64d..2230c79f 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -21,6 +21,7 @@ using GeophysicalFlows, CUDA, Random, Printf, CairoMakie parsevalsum = FourierFlows.parsevalsum +record = CairoMakie.record # ## Choosing a device: CPU or GPU From e3bf340bec380e4cf25b61dbb2df3693bd7ecce8 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 16:19:04 +1000 Subject: [PATCH 32/44] add Lx, Ly shortcuts --- examples/twodnavierstokes_stochasticforcing_budgets.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 66920a4b..2209ea58 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -97,7 +97,8 @@ 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 From e0450409093662f14c9e9bb09903948f42d8cb60 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 17:04:11 +1000 Subject: [PATCH 33/44] literated commented code --- examples/multilayerqg_2layer.jl | 6 +++--- examples/surfaceqg_decaying.jl | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 30dfe713..7a1677c6 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -146,13 +146,13 @@ q₂ = Observable(Array(vars.q[:, :, 2])) ψ₂ = Observable(Array(vars.ψ[:, :, 2])) function compute_levels(maxf, nlevels=8) - # -max(|f|):...:max(|f|) + ## -max(|f|):...:max(|f|) levelsf = @lift collect(range(-$maxf, stop = $maxf, length=nlevels)) - # only positive + ## only positive levelsf⁺ = @lift collect(range($maxf/(nlevels-1), stop = $maxf, length=Int(nlevels/2))) - # only negative + ## only negative levelsf⁻ = @lift collect(range(-$maxf, stop = -$maxf/(nlevels-1), length=Int(nlevels/2))) return levelsf, levelsf⁺, levelsf⁻ diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index f463eaba..9355c064 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -207,14 +207,14 @@ end fig = Figure(resolution = (800, 380)) -kwargs_axis = (xlabel = "x", +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) * ")", kwargs_axis...) -axu = Axis(fig[1, 2]; title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", kwargs_axis...) -axv = Axis(fig[1, 3]; title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", kwargs_axis...) +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, vars.b; colormap = :deep, colorrange = (0, 1)) From 0a8962c9de6f5cc51f1ee0df44a16131d1be9a5d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 17:25:25 +1000 Subject: [PATCH 34/44] gpu friendly --- examples/barotropicqgql_betaforced.jl | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 65b0e64f..71b8ddb2 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -112,6 +112,7 @@ 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) fig = Figure() @@ -123,7 +124,7 @@ ax = Axis(fig[1, 1], title = "a forcing realization", limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) -heatmap!(ax, x, y, irfft(vars.Fh, grid.nx); +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); colormap = :balance, colorrange = (-8, 8)) fig @@ -227,11 +228,11 @@ axZ = Axis(fig[2, 3], limits = ((-0.1, 4.1), (0, 5))) ζ̄, ζ′= prob.vars.Zeta, prob.vars.zeta -ζ = Observable(@. ζ̄ + ζ′) +ζ = Observable(Array(@. ζ̄ + ζ′)) ψ̄, ψ′= prob.vars.Psi, prob.vars.psi -ψ = Observable(@. ψ̄ + ψ′) -ζ̄ₘ = Observable(vec(mean(ζ̄, dims=1))) -ūₘ = Observable(vec(mean(prob.vars.U, dims=1))) +ψ = 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]) @@ -274,10 +275,10 @@ record(fig, "barotropicqgql_betaforced.mp4", frames, framerate = 18) do j println(log) end - ζ[] = @. ζ̄ + ζ′ - ψ[] = @. ψ̄ + ψ′ - ζ̄ₘ[] = vec(mean(ζ̄, dims=1)) - ūₘ[] = vec(mean(prob.vars.U, dims=1)) + ζ[] = Array(@. ζ̄ + ζ′) + ψ[] = Array(@. ψ̄ + ψ′) + ζ̄ₘ[] = Array(vec(mean(ζ̄, dims=1))) + ūₘ[] = Array(vec(mean(prob.vars.U, dims=1))) μt.val = μ * E.t[1:E.i] energy[] = E.data[1:E.i] From 5ab97bcb24a7767a1753d31378e304f440a5b18d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 17:41:34 +1000 Subject: [PATCH 35/44] ensure gpu friendly --- examples/barotropicqgql_betaforced.jl | 8 ++++---- examples/multilayerqg_2layer.jl | 8 ++++---- examples/singlelayerqg_betadecay.jl | 10 +++++----- examples/singlelayerqg_betaforced.jl | 6 +++--- ...layerqg_decaying_barotropic_equivalentbarotropic.jl | 4 ++-- examples/singlelayerqg_decaying_topography.jl | 10 +++++----- examples/surfaceqg_decaying.jl | 10 +++++----- examples/twodnavierstokes_decaying.jl | 4 ++-- examples/twodnavierstokes_stochasticforcing.jl | 6 +++--- examples/twodnavierstokes_stochasticforcing_budgets.jl | 2 +- 10 files changed, 34 insertions(+), 34 deletions(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 71b8ddb2..96332e8d 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -275,10 +275,10 @@ record(fig, "barotropicqgql_betaforced.mp4", frames, framerate = 18) do j println(log) end - ζ[] = Array(@. ζ̄ + ζ′) - ψ[] = Array(@. ψ̄ + ψ′) - ζ̄ₘ[] = Array(vec(mean(ζ̄, dims=1))) - ūₘ[] = Array(vec(mean(prob.vars.U, dims=1))) + ζ[] = @. ζ̄ + ζ′ + ψ[] = @. ψ̄ + ψ′ + ζ̄ₘ[] = vec(mean(ζ̄, dims=1)) + ūₘ[] = vec(mean(prob.vars.U, dims=1)) μt.val = μ * E.t[1:E.i] energy[] = E.data[1:E.i] diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 30dfe713..b2a5aaa7 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -241,10 +241,10 @@ record(fig, "multilayerqg_2layer.mp4", frames, framerate = 18) do j println(log) end - q₁[] = Array(vars.q[:, :, 1]) - ψ₁[] = Array(vars.ψ[:, :, 1]) - q₂[] = Array(vars.q[:, :, 2]) - ψ₂[] = Array(vars.ψ[:, :, 2]) + q₁[] = vars.q[:, :, 1] + ψ₁[] = vars.ψ[:, :, 1] + q₂[] = vars.q[:, :, 2] + ψ₂[] = vars.ψ[:, :, 2] maxψ₁[] = maximum(abs, vars.ψ[:, :, 1]) maxψ₂[] = maximum(abs, vars.ψ[:, :, 2]) diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 5742fa3a..e6b28968 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -23,7 +23,7 @@ using Statistics: mean # ## Choosing a device: CPU or GPU -dev = CPU() # Device (CPU/GPU) +dev = GPU() # Device (CPU/GPU) nothing # hide @@ -171,10 +171,10 @@ axū = Axis(fig[2, 2], aspect = 1, limits = ((-0.5, 0.5), (-Ly/2, Ly/2))) -q = Observable(vars.q) -ψ = Observable(vars.ψ) -q̄ₘ = Observable(vec(mean(vars.q, dims=1))) -ūₘ = Observable(vec(mean(vars.u, dims=1))) +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)) diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index aa89faff..395332da 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -159,7 +159,7 @@ if !isdir(plotpath); mkdir(plotpath); end nothing # hide # and then create Output. -get_sol(prob) = prob.sol # extracts the Fourier-transformed solution +get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution function get_u(prob) vars, grid, sol = prob.vars, prob.grid, prob.sol @@ -170,7 +170,7 @@ function get_u(prob) ldiv!(vars.u, grid.rfftplan, -im * grid.l .* vars.ψh) - return vars.u + return Array(vars.u) end output = Output(prob, filename, (:qh, get_sol), (:u, get_u)) @@ -237,7 +237,7 @@ Lx, Ly = file["grid/Lx"], file["grid/Ly"] j = Observable(1) q = @lift irfft(qh[$j], nx) -ψ = @lift irfft(- grid.invKrsq .* 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)) diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index c8bba245..d54907d7 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -121,8 +121,8 @@ title_eqbqg = @lift "equivalent barotropic; deformation radius: " * @sprintf("%. ax1 = Axis(fig[1, 1]; title = title_bqg, axis_kwargs...) ax2 = Axis(fig[1, 2]; title = title_eqbqg, axis_kwargs...) -ζ_bqg = Observable(relativevorticity(prob_bqg)) -ζ_eqbqg = Observable(relativevorticity(prob_eqbqg)) +ζ_bqg = Observable(Array(relativevorticity(prob_bqg))) +ζ_eqbqg = Observable(Array(relativevorticity(prob_eqbqg))) heatmap!(ax1, x, y, ζ_bqg; colormap = :balance, colorrange = (-40, 40)) diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 241f02cd..f79f48c3 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -70,6 +70,8 @@ nothing # hide # plotted with `Array()` to make sure it is brought back on the CPU when the variable lives # on the GPU. +η = Array(params.eta) + fig = Figure() ax = Axis(fig[1, 1]; xlabel = "x", @@ -77,7 +79,7 @@ ax = Axis(fig[1, 1]; title = "topographic PV η=f₀h/H", limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) -contourf!(ax, x, y, params.eta; +contourf!(ax, x, y, η; levels = collect(-3:0.4:3), colormap = :balance, colorrange = (-3, 3)) fig @@ -105,8 +107,8 @@ nothing # hide # Let's plot the initial vorticity and streamfunction. -q = Observable(vars.q) -ψ = Observable(vars.ψ) +q = Observable(Array(vars.q)) +ψ = Observable(Array(vars.ψ)) fig = Figure(resolution=(800, 380)) @@ -170,8 +172,6 @@ nothing # hide # We modify the figure with the initial state slightly. We add the topography contours and # also the time. -η = prob.params.eta - contour!(axq, x, y, η; levels = collect(0.5:0.5:3), linewidth = 2, color = (:black, 0.5)) diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index f463eaba..05e60d46 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -139,7 +139,7 @@ nothing # hide # We define a function that plots the buoyancy field and the time evolution of kinetic energy # and buoyancy variance. -b = Observable(vars.b) +b = Observable(Array(vars.b)) ke = Observable([Point2f(KE.t[1], KE.data[1])]) b² = Observable([Point2f(B.t[1], B.data[1])]) @@ -189,7 +189,7 @@ record(fig, "sqg_ellipticalvortex.mp4", 0:round(Int, nsteps/nsubs), framerate = println(log2) end - b[] = Array(vars.b) + 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])) @@ -216,17 +216,17 @@ axb = Axis(fig[1, 1]; title = "bₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", axu = Axis(fig[1, 2]; title = "uₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", kwargs_axis...) axv = Axis(fig[1, 3]; title = "vₛ(x, y, t=" * @sprintf("%.2f", clock.t) * ")", kwargs_axis...) -hb = heatmap!(axb, x, y, vars.b; +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, vars.u; +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, vars.v; +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) diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index b5567be2..fdd5396d 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -122,7 +122,7 @@ nothing # hide # We initialize a plot with the vorticity field and the time-series of # energy and enstrophy diagnostics. -ζ = Observable(vars.ζ) +ζ = Observable(Array(vars.ζ)) title_ζ = Observable("vorticity, t=" * @sprintf("%.2f", clock.t)) energy = Observable(Point2f[(E.t[1], E.data[1] / E.data[1])]) @@ -166,7 +166,7 @@ record(fig, "twodturb.mp4", 0:Int(nsteps/nsubs), framerate = 18) do j println(log) end - ζ[] = Array(vars.ζ) + ζ[] = 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])) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 2230c79f..accc041d 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -116,7 +116,7 @@ ax = Axis(fig[1, 1], title = "a forcing realization", limits = ((-L/2, L/2), (-L/2, L/2))) -heatmap!(ax, x, y, irfft(vars.Fh, grid.nx); +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); colormap = :balance, colorrange = (-200, 200)) fig @@ -143,7 +143,7 @@ nothing # hide # energy and enstrophy diagnostics. To plot energy and enstrophy on the # same axes we scale enstrophy with ``k_f^2``. -ζ = Observable(vars.ζ) +ζ = Observable(Array(vars.ζ)) title_ζ = Observable("vorticity, μ t=" * @sprintf("%.2f", μ * clock.t)) energy = Observable(Point2f[(μ * E.t[1], E.data[1])]) @@ -187,7 +187,7 @@ record(fig, "twodturb_forced.mp4", 0:round(Int, nsteps / nsubs), framerate = 18) println(log) end - ζ[] = Array(vars.ζ) + ζ[] = 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)) diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index 2209ea58..eefa0814 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -119,7 +119,7 @@ ax = Axis(fig[1, 1], title = "a forcing realization", limits = ((-Lx/2, Lx/2), (-Ly/2, Ly/2))) -heatmap!(ax, x, y, irfft(vars.Fh, grid.nx); +heatmap!(ax, x, y, Array(irfft(vars.Fh, grid.nx)); colormap = :balance, colorrange = (-200, 200)) fig From 865f66af50ab6c6c3a466bdabf3569557fce5ba3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 17:43:00 +1000 Subject: [PATCH 36/44] add example refs --- docs/src/modules/barotropicqgql.md | 2 +- docs/src/modules/multilayerqg.md | 2 +- docs/src/modules/singlelayerqg.md | 8 ++++---- docs/src/modules/surfaceqg.md | 2 +- docs/src/modules/twodnavierstokes.md | 6 +++--- 5 files changed, 10 insertions(+), 10 deletions(-) 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..d463266f 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, producing an animation of the evolution of the surface buoyancy. > 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. From 78b7734747d8020da62d0aa3f7bdf4f299ec9a0d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 17:45:04 +1000 Subject: [PATCH 37/44] drop GeophysicalFlows from docs/Project --- docs/Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 97319fb9..e0a67d1e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,7 +2,6 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -GeophysicalFlows = "44ee3b1c-bc02-53fa-8355-8e347616e15e" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" From 360feeb3081565e5519f15dfd01962084753c667 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 17:49:40 +1000 Subject: [PATCH 38/44] back to CPU --- examples/singlelayerqg_betadecay.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index e6b28968..7aebca3d 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -23,7 +23,7 @@ using Statistics: mean # ## Choosing a device: CPU or GPU -dev = GPU() # Device (CPU/GPU) +dev = CPU() # Device (CPU/GPU) nothing # hide From 729a098b23af8aae5e049ace0fbae861e0211837 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 18:50:14 +1000 Subject: [PATCH 39/44] some cosmetics --- examples/barotropicqgql_betaforced.jl | 1 + examples/multilayerqg_2layer.jl | 1 + examples/singlelayerqg_betadecay.jl | 1 + examples/singlelayerqg_betaforced.jl | 1 + ...ecaying_barotropic_equivalentbarotropic.jl | 5 ++++- examples/singlelayerqg_decaying_topography.jl | 1 + examples/surfaceqg_decaying.jl | 5 +++-- examples/twodnavierstokes_decaying.jl | 19 ++++++++----------- .../twodnavierstokes_stochasticforcing.jl | 1 + 9 files changed, 21 insertions(+), 14 deletions(-) diff --git a/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 96332e8d..3d2072c2 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -289,6 +289,7 @@ record(fig, "barotropicqgql_betaforced.mp4", frames, framerate = 18) do j stepforward!(prob, diags, nsubs) BarotropicQGQL.updatevars!(prob) end +nothing # hide # ![](barotropicqgql_betaforced.mp4) diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index 05f5deaa..be3cb1d2 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -258,6 +258,7 @@ record(fig, "multilayerqg_2layer.mp4", frames, framerate = 18) do j stepforward!(prob, diags, nsubs) MultiLayerQG.updatevars!(prob) end +nothing # hide # ![](multilayerqg_2layer.mp4) diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 7aebca3d..751d97e9 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -223,6 +223,7 @@ record(fig, "singlelayerqg_betadecay.mp4", frames, framerate = 8) do j stepforward!(prob, diags, nsubs) SingleLayerQG.updatevars!(prob) end +nothing # hide # ![](singlelayerqg_betadecay.mp4) diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 395332da..bff838ed 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -312,5 +312,6 @@ record(fig, "singlelayerqg_betaforced.mp4", frames, framerate = 18) do 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 d54907d7..6cee3ef4 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -92,7 +92,9 @@ SingleLayerQG.set_q!(prob_eqbqg, q₀_eqbqg) nothing # hide -# Let's plot the initial vorticity field for each problem. +# 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 @@ -162,5 +164,6 @@ record(fig, "singlelayerqg_barotropic_equivalentbarotropic.mp4", 0:Int(nsteps/ns ζ_bqg[] = relativevorticity(prob_bqg) ζ_eqbqg[] = relativevorticity(prob_eqbqg) end +nothing # hide # ![](singlelayerqg_barotropic_equivalentbarotropic.mp4) diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index f79f48c3..07edb2d8 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -209,5 +209,6 @@ record(fig, "singlelayerqg_decaying_topography.mp4", 0:round(Int, nsteps/nsubs), stepforward!(prob, diags, nsubs) SingleLayerQG.updatevars!(prob) end +nothing # hide # ![](singlelayerqg_decaying_topography.mp4) diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index a52bd6a6..3a27bdc7 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -199,6 +199,7 @@ record(fig, "sqg_ellipticalvortex.mp4", 0:round(Int, nsteps/nsubs), framerate = stepforward!(prob, diags, nsubs) SurfaceQG.updatevars!(prob) end +nothing # hide # ![](sqg_ellipticalvortex.mp4) @@ -236,7 +237,7 @@ fig # ## Save -# Last we can save the output by calling +# Last, we can save the output by calling # ```julia -# saveoutput(out)` +# saveoutput(out) # ``` diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index fdd5396d..a4fa1c52 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -176,31 +176,28 @@ record(fig, "twodturb.mp4", 0:Int(nsteps/nsubs), framerate = 18) do j stepforward!(prob, diags, nsubs) TwoDNavierStokes.updatevars!(prob) end +nothing # hide # ![](twodturb.mp4) -# Last we can save the output by calling -# ```julia -# saveoutput(out)` -# ``` - - # ## Radial energy spectrum # 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. lines(kr, vec(abs.(Ehr)); linewidth = 2, - axis = (xlabel = "kᵣ", - ylabel = "∫ |Ê| kᵣ dk_θ", + axis = (xlabel = L"k_r", + ylabel = L"\int |\hat{E}| k_r \mathrm{d}k_\theta", xscale = log10, yscale = log10, title = "Radial energy spectrum", diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index accc041d..cdae9ba2 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -197,5 +197,6 @@ record(fig, "twodturb_forced.mp4", 0:round(Int, nsteps / nsubs), framerate = 18) stepforward!(prob, diags, nsubs) TwoDNavierStokes.updatevars!(prob) end +nothing # hide # ![](twodturb_forced.mp4) From 76aa3855e4be194f244c642cee90f9f3ec83c366 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 21:18:33 +1000 Subject: [PATCH 40/44] require FourierFlows v0.9.4 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index e6ccc45f..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.3" +FourierFlows = "^0.9.4" JLD2 = "^0.1, ^0.2, ^0.3, ^0.4" Reexport = "^0.2, ^1" SpecialFunctions = "^0.10, ^1, 2" From aa029ca7cd2c33b2b0959f38b0fae6f0ce4bb468 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 21:19:28 +1000 Subject: [PATCH 41/44] add suggestions by @jbisits --- docs/src/modules/surfaceqg.md | 2 +- examples/barotropicqgql_betaforced.jl | 2 +- examples/multilayerqg_2layer.jl | 2 +- examples/singlelayerqg_betadecay.jl | 4 ++-- examples/singlelayerqg_betaforced.jl | 2 +- ...layerqg_decaying_barotropic_equivalentbarotropic.jl | 2 +- examples/singlelayerqg_decaying_topography.jl | 10 +++++----- examples/surfaceqg_decaying.jl | 4 ++-- examples/twodnavierstokes_decaying.jl | 2 +- examples/twodnavierstokes_stochasticforcing.jl | 4 ++-- examples/twodnavierstokes_stochasticforcing_budgets.jl | 6 +++--- 11 files changed, 20 insertions(+), 20 deletions(-) diff --git a/docs/src/modules/surfaceqg.md b/docs/src/modules/surfaceqg.md index d463266f..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`](@ref surfaceqg_decaying_example): Simulate 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/examples/barotropicqgql_betaforced.jl b/examples/barotropicqgql_betaforced.jl index 3d2072c2..23765e51 100644 --- a/examples/barotropicqgql_betaforced.jl +++ b/examples/barotropicqgql_betaforced.jl @@ -258,7 +258,7 @@ nothing # hide # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() diff --git a/examples/multilayerqg_2layer.jl b/examples/multilayerqg_2layer.jl index be3cb1d2..68caa332 100644 --- a/examples/multilayerqg_2layer.jl +++ b/examples/multilayerqg_2layer.jl @@ -16,7 +16,7 @@ # ``` # ## 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, CairoMakie, Printf diff --git a/examples/singlelayerqg_betadecay.jl b/examples/singlelayerqg_betadecay.jl index 751d97e9..3d275ce0 100644 --- a/examples/singlelayerqg_betadecay.jl +++ b/examples/singlelayerqg_betadecay.jl @@ -14,7 +14,7 @@ # ``` # ## 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, CairoMakie, Printf, Random @@ -197,7 +197,7 @@ fig # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index bff838ed..86d66373 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -15,7 +15,7 @@ # ``` # ## 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, JLD2, CairoMakie, Random, Printf diff --git a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl index 6cee3ef4..b2775673 100644 --- a/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl +++ b/examples/singlelayerqg_decaying_barotropic_equivalentbarotropic.jl @@ -15,7 +15,7 @@ # ``` # ## 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, CairoMakie diff --git a/examples/singlelayerqg_decaying_topography.jl b/examples/singlelayerqg_decaying_topography.jl index 07edb2d8..0a320ba8 100644 --- a/examples/singlelayerqg_decaying_topography.jl +++ b/examples/singlelayerqg_decaying_topography.jl @@ -10,11 +10,11 @@ # ```julia # using Pkg -# pkg"add GeophysicalFlows, CairoMakie, 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, CairoMakie, Printf, Random @@ -169,8 +169,8 @@ nothing # hide # ## Visualizing the simulation -# We modify the figure with the initial state slightly. We add the topography contours and -# also the time. +# 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)) @@ -186,7 +186,7 @@ nothing # hide # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() diff --git a/examples/surfaceqg_decaying.jl b/examples/surfaceqg_decaying.jl index 3a27bdc7..5e95accb 100644 --- a/examples/surfaceqg_decaying.jl +++ b/examples/surfaceqg_decaying.jl @@ -19,7 +19,7 @@ # ``` # ## 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, CairoMakie, Printf, Random @@ -237,7 +237,7 @@ fig # ## 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) # ``` diff --git a/examples/twodnavierstokes_decaying.jl b/examples/twodnavierstokes_decaying.jl index a4fa1c52..2ad6efea 100644 --- a/examples/twodnavierstokes_decaying.jl +++ b/examples/twodnavierstokes_decaying.jl @@ -14,7 +14,7 @@ # ``` # ## 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, CairoMakie diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index cdae9ba2..b0fc5abd 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -16,7 +16,7 @@ # ``` # ## 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 @@ -174,7 +174,7 @@ fig # ## Time-stepping the `Problem` forward -# Finally, we time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() diff --git a/examples/twodnavierstokes_stochasticforcing_budgets.jl b/examples/twodnavierstokes_stochasticforcing_budgets.jl index eefa0814..7707f475 100644 --- a/examples/twodnavierstokes_stochasticforcing_budgets.jl +++ b/examples/twodnavierstokes_stochasticforcing_budgets.jl @@ -18,7 +18,7 @@ # ``` # ## 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 @@ -147,7 +147,7 @@ nothing # hide # ## Time-stepping the `Problem` forward -# We time-step the `Problem` forward in time. +# We step the `Problem` forward in time. startwalltime = time() for i = 1:ns @@ -181,7 +181,7 @@ heatmap!(ax, x, y, Array(vars.ζ); fig -# And finaly, we plot the evolution of the energy and enstrophy diagnostics and all terms +# 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``. From 36ffd97f77e98589963879fcc6c8e8cb8a82748e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 21:26:59 +1000 Subject: [PATCH 42/44] minor change to trigger CI --- README.md | 1 - 1 file changed, 1 deletion(-) diff --git a/README.md b/README.md index b4248694..90c42090 100644 --- a/README.md +++ b/README.md @@ -143,7 +143,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 From 78295e470c298085238990cf0b7ccf795943067d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 21:28:16 +1000 Subject: [PATCH 43/44] minor phrasing change --- examples/singlelayerqg_betaforced.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/singlelayerqg_betaforced.jl b/examples/singlelayerqg_betaforced.jl index 86d66373..bbceac9d 100644 --- a/examples/singlelayerqg_betaforced.jl +++ b/examples/singlelayerqg_betaforced.jl @@ -212,8 +212,8 @@ savediagnostic(Z, "enstrophy", output.path) # ## Load saved output and visualize -# We now have output from our simulation saved in `singlelayerqg_forcedbeta.jld2`. We can -# now load the JLD2 output and create a time series for fields that interest us. +# 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) From e0a1975d2e9df85a62702e6374c6873e241be5c4 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 30 Aug 2022 21:47:47 +1000 Subject: [PATCH 44/44] some rephrasing --- README.md | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 90c42090..e8d4ced5 100644 --- a/README.md +++ b/README.md @@ -62,14 +62,17 @@ Some animations created with GeophysicalFlows.jl are [online @ youtube]. ## Modules -* [`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 and with the ability - 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. +* [`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