From 245dd7c0254c98d4b9a421f4fdc73867fe9735a7 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Mon, 20 Jul 2020 15:59:37 -0700 Subject: [PATCH 01/19] Add Mat file writer --- .gitignore | 1 + examples/2D_Paper/2D_Paper_Inverse_Cascade.jl | 311 ++++++++++++++++++ 2 files changed, 312 insertions(+) create mode 100644 examples/2D_Paper/2D_Paper_Inverse_Cascade.jl diff --git a/.gitignore b/.gitignore index 18027391..5cd30ca1 100644 --- a/.gitignore +++ b/.gitignore @@ -38,3 +38,4 @@ coverage/ # MacOS stuff *.DS_Store +*.mat diff --git a/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl b/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl new file mode 100644 index 00000000..bc634ecd --- /dev/null +++ b/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl @@ -0,0 +1,311 @@ +# # 2D forced-dissipative turbulence +# +#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). +# +# A simulation of forced-dissipative two-dimensional turbulence. We solve the +# two-dimensional vorticity equation with stochastic excitation and dissipation in +# the form of linear drag and hyperviscosity. As a demonstration, we compute how +# each of the forcing and dissipation terms contribute to the energy budget. + +using FourierFlows, Printf, Random, Plots + +using FourierFlows: parsevalsum +using Random: seed! +using FFTW: irfft + +import GeophysicalFlows.TwoDNavierStokes +import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy, dissipation, work, drag +import GeophysicalFlows: peakedisotropicspectrum + + +# ## Choosing a device: CPU or GPU + +dev = CPU() # Device (CPU/GPU) +nothing # hide + + +# ## Numerical, domain, and simulation parameters +# +# First, we pick some numerical and physical parameters for our model. + + n, L = 256, 2π # grid resolution and domain length + ν, nν = 2e-7, 2 # hyperviscosity coefficient and order + μ, nμ = 1e-1, 0 # linear drag coefficient +dt, tf = 0.005, 0.2/μ # timestep and final time + nt = round(Int, tf/dt) # total timesteps + ns = 4 # how many intermediate times we want to plot +nothing # hide + + +# ## Forcing + +# We force the vorticity equation with stochastic excitation that is delta-correlated +# in time and while spatially homogeneously and isotropically correlated. The forcing +# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and +# width $\delta k_f$, and it injects energy per unit area and per unit time equal +# to $\varepsilon$. + +forcing_wavenumber = 14.0 # the central forcing wavenumber for a spectrum that is a ring in wavenumber space +forcing_bandwidth = 1.5 # the width of the forcing spectrum +ε = 0.001 # energy input rate by the forcing + +gr = TwoDGrid(dev, n, L) +x, y = gr.x, gr.y + +forcing_spectrum = @. exp(-(sqrt(gr.Krsq)-forcing_wavenumber)^2/(2*forcing_bandwidth^2)) +forcing_spectrum[ gr.Krsq .< (2π/L*2)^2 ] .= 0 # make sure that focing has no power for low wavenumbers +forcing_spectrum[ gr.Krsq .> (2π/L*20)^2 ] .= 0 # make sure that focing has no power for high wavenumbers +ε0 = parsevalsum(forcing_spectrum.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly) +forcing_spectrum .= ε/ε0 * forcing_spectrum # normalize forcing to inject energy ε + +seed!(1234) +nothing # hide + +# Next we construct function `calcF!` that computes a forcing realization every timestep +function calcF!(Fh, sol, t, cl, v, p, g) + ξ = ArrayType(dev)(exp.(2π*im*rand(eltype(gr), size(sol)))/sqrt(cl.dt)) + ξ[1, 1] = 0 + @. Fh = ξ*sqrt(forcing_spectrum) + nothing +end +nothing # hide + + +# ## Problem setup +# We initialize a `Problem` by providing a set of keyword arguments. The +# `stepper` keyword defines the time-stepper to be used. +prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4", + calcF=calcF!, stochastic=true) +nothing # hide + +# Define some shortcuts for convenience. +sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid +nothing # hide + + +# First let's see how a forcing realization looks like. +calcF!(v.Fh, sol, 0.0, cl, v, p, g) + +heatmap(x, y, irfft(v.Fh, g.nx), + aspectratio = 1, + c = :balance, + clim = (-200, 200), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "x", + ylabel = "y", + title = "a forcing realization", + framestyle = :box) + + +# ## Setting initial conditions + +# Our initial condition is simply fluid at rest. +TwoDNavierStokes.set_zeta!(prob, zeros(g.nx, g.ny)) + + +# ## Diagnostics + +# Create Diagnostics; the diagnostics are aimed to probe the energy budget. +E = Diagnostic(energy, prob, nsteps=nt) # energy +R = Diagnostic(drag, prob, nsteps=nt) # dissipation by drag +D = Diagnostic(dissipation, prob, nsteps=nt) # dissipation by hyperviscosity +W = Diagnostic(work, prob, nsteps=nt) # work input by forcing +diags = [E, D, W, R] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. +nothing # hide + + +# ## Visualizing the simulation + +# We define a function that plots the vorticity field and the evolution of +# the diagnostics: energy and all terms involved in the energy budget. Last +# we confirm whether the energy budget is accurate, i.e., $\mathrm{d}E/\mathrm{d}t = W - R - D$. + +function computetendencies_and_makeplot(prob, diags) + TwoDNavierStokes.updatevars!(prob) + E, D, W, R = diags + + clocktime = round(μ*cl.t, digits=2) + + i₀ = 1 + dEdt_numerical = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt #numerical first-order approximation of energy tendency + ii = (i₀):E.i-1 + ii2 = (i₀+1):E.i + + t = E.t[ii] + dEdt_computed = W[ii2] - D[ii] - R[ii] # Stratonovich interpretation + + residual = dEdt_computed - dEdt_numerical + + l = @layout grid(2, 2) + + p1 = heatmap(x, y, v.zeta, + aspectratio = 1, + legend = false, + c = :viridis, + clim = (-25, 25), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "μt", + ylabel = "y", + title = "∇²ψ(x, y, t="*@sprintf("%.2f", cl.t)*")", + framestyle = :box) + + p2 = plot(μ*t, [W[ii2] ε.+0*t -D[ii] -R[ii]], + label = ["work, W" "ensemble mean work, " "dissipation, D" "drag, D=-2μE"], + linestyle = [:solid :dash :solid :solid], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "energy sources and sinks") + + p3 = plot(μ*t, [dEdt_computed[ii], dEdt_numerical], + label = ["computed W-D" "numerical dE/dt"], + linestyle = [:solid :dashdotdot], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "dE/dt") + + p4 = plot(μ*t, residual, + label = "residual dE/dt = computed - numerical", + linewidth = 2, + alpha = 0.7, + xlabel = "μt") + + p = plot(p1, p2, p3, p4, layout=l, size = (900, 800)) + return p +end +nothing # hide + + + + +# ## Time-stepping the `Problem` forward + +# Finally, we time-step the `Problem` forward in time. + +startwalltime = time() +for i = 1:ns + stepforward!(prob, diags, round(Int, nt/ns)) + TwoDNavierStokes.updatevars!(prob) + cfl = cl.dt*maximum([maximum(v.u)/g.dx, maximum(v.v)/g.dy]) + + log = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min", cl.step, cl.t, + cfl, (time()-startwalltime)/60) + + println(log) +end + + +# ## Plot +# And now let's see what we got. We plot the output. + +p = computetendencies_and_makeplot(prob, diags) + +# ## Output + +# We choose folder for outputing `.jld2` files and snapshots (`.png` files). +filepath = "." +plotpath = "/Users/brodiepearson/GitHub/GeophysicalFlows.jl/examples/2D_Paper/Plots/Inverse_Cascade" +plotname = "snapshots" +filename = joinpath(filepath, "Inverse_Cascade.jld2") +nothing # hide + +# Do some basic file management +if isfile(filename); rm(filename); end +if !isdir(plotpath); mkdir(plotpath); end +nothing # hide + +# And then create Output +get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution +get_u(prob) = Array(irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)) +out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) +saveproblem(out) +nothing # hide + +# ## Visualizing the simulation + +# We initialize a plot with the vorticity field and the time-series of +# energy and enstrophy diagnostics. + +p1 = heatmap(x, y, v.zeta, + aspectratio = 1, + c = :balance, + clim = (-40, 40), + xlims = (-L/2, L/2), + ylims = (-L/2, L/2), + xticks = -3:3, + yticks = -3:3, + xlabel = "x", + ylabel = "y", + title = "vorticity, t="*@sprintf("%.2f", cl.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 grid(1, 2) +p = plot(p1, p2, layout = l, size = (900, 400)) + + +# ## Time-stepping the `Problem` forward + +# We time-step the `Problem` forward in time. + +startwalltime = time() + +anim = @animate for j = 0:Int(nsteps/nsubs) + + log = @sprintf("step: %04d, t: %d, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min", + cl.step, cl.t, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60) + + if j%(1000/nsubs)==0; println(log) end + + p[1][1][:z] = vs.zeta + p[1][:title] = "vorticity, t="*@sprintf("%.2f", cl.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]) + + stepforward!(prob, diags, nsubs) + TwoDNavierStokes.updatevars!(prob) + +end + +mp4(anim, "twodturb.mp4", fps=18) + + +# Last we save the output. +saveoutput(out) + + +# ## Radial energy spectrum + +# After the simulation is done we plot the radial energy spectrum to illustrate +# how `FourierFlows.radialspectrum` can be used, + +E = @. 0.5*(vs.u^2 + vs.v^2) # energy density +Eh = rfft(E) # Fourier transform of energy density +kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1) # compute radial specturm of `Eh` +nothing # hide + +# and we plot it. +plot(kr, abs.(Ehr), + linewidth = 2, + alpha = 0.7, + xlabel = "kᵣ", ylabel = "∫ |Ê| kᵣ dk_θ", + xlims = (5e-1, gr.nx), + xscale = :log10, yscale = :log10, + title = "Radial energy spectrum", + legend = false) From c716b58dfb50a4f295a613f781610fa418ce75bc Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Mon, 20 Jul 2020 16:32:02 -0700 Subject: [PATCH 02/19] Update 2D_Paper_Inverse_Cascade.jl --- examples/2D_Paper/2D_Paper_Inverse_Cascade.jl | 33 ++++++++++++++----- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl b/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl index bc634ecd..52a77944 100644 --- a/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl +++ b/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl @@ -12,7 +12,9 @@ using FourierFlows, Printf, Random, Plots using FourierFlows: parsevalsum using Random: seed! -using FFTW: irfft +using FFTW: rfft, irfft + +using MAT import GeophysicalFlows.TwoDNavierStokes import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy, dissipation, work, drag @@ -34,7 +36,8 @@ nothing # hide μ, nμ = 1e-1, 0 # linear drag coefficient dt, tf = 0.005, 0.2/μ # timestep and final time nt = round(Int, tf/dt) # total timesteps - ns = 4 # how many intermediate times we want to plot + ns = 4 # how many intermediate times we want to plot + nsubs = round(Int, nt/ns) nothing # hide @@ -113,8 +116,9 @@ TwoDNavierStokes.set_zeta!(prob, zeros(g.nx, g.ny)) E = Diagnostic(energy, prob, nsteps=nt) # energy R = Diagnostic(drag, prob, nsteps=nt) # dissipation by drag D = Diagnostic(dissipation, prob, nsteps=nt) # dissipation by hyperviscosity -W = Diagnostic(work, prob, nsteps=nt) # work input by forcing -diags = [E, D, W, R] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. +W = Diagnostic(work, prob, nsteps=nt) # work input by +Z = Diagnostic(enstrophy, prob; nsteps=nt) +diags = [E, D, W, R, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -126,7 +130,7 @@ nothing # hide function computetendencies_and_makeplot(prob, diags) TwoDNavierStokes.updatevars!(prob) - E, D, W, R = diags + E, D, W, R, Z = diags clocktime = round(μ*cl.t, digits=2) @@ -215,6 +219,7 @@ filepath = "." plotpath = "/Users/brodiepearson/GitHub/GeophysicalFlows.jl/examples/2D_Paper/Plots/Inverse_Cascade" plotname = "snapshots" filename = joinpath(filepath, "Inverse_Cascade.jld2") +matfilename = joinpath(filepath, "Inverse_Cascade.mat") nothing # hide # Do some basic file management @@ -227,6 +232,18 @@ get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution get_u(prob) = Array(irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)) out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) saveproblem(out) +matfile = matopen(matfilename) +u_temp = (:u, get_u) +write(matfile, "zeta", v.zeta) +write(matfile, "u", v.u) +write(matfile, "v", v.v) +write(matfile, "nx", gr.nx) +write(matfile, "ny", gr.ny) +write(matfile, "dx", gr.dx) +write(matfile, "dy", gr.dy) +write(matfile, "Lx", gr.Lx) +write(matfile, "Ly", gr.Ly) +close(matfile) nothing # hide # ## Visualizing the simulation @@ -266,14 +283,14 @@ p = plot(p1, p2, layout = l, size = (900, 400)) startwalltime = time() -anim = @animate for j = 0:Int(nsteps/nsubs) +anim = @animate for j = 0:Int(nt/nsubs) log = @sprintf("step: %04d, t: %d, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min", cl.step, cl.t, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60) if j%(1000/nsubs)==0; println(log) end - p[1][1][:z] = vs.zeta + p[1][1][:z] = v.zeta p[1][:title] = "vorticity, t="*@sprintf("%.2f", cl.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]) @@ -295,7 +312,7 @@ saveoutput(out) # After the simulation is done we plot the radial energy spectrum to illustrate # how `FourierFlows.radialspectrum` can be used, -E = @. 0.5*(vs.u^2 + vs.v^2) # energy density +E = @. 0.5*(v.u^2 + v.v^2) # energy density Eh = rfft(E) # Fourier transform of energy density kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1) # compute radial specturm of `Eh` nothing # hide From 029f1d363262936fe560364cc1891f9057020ebe Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Wed, 29 Jul 2020 10:10:09 -0700 Subject: [PATCH 03/19] Update twodnavierstokes.jl --- src/twodnavierstokes.jl | 25 +++++++++++++++++++------ 1 file changed, 19 insertions(+), 6 deletions(-) diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl index 1796866d..181f8aa0 100644 --- a/src/twodnavierstokes.jl +++ b/src/twodnavierstokes.jl @@ -9,7 +9,8 @@ export enstrophy, dissipation, work, - drag + drag, + drag_ens using Reexport @@ -44,13 +45,13 @@ function Problem(dev::Device=CPU(); T = Float64) grid = TwoDGrid(dev, nx, Lx, ny, Ly; T=T) - + params = Params{T}(ν, nν, μ, nμ, calcF) - + vars = calcF == nothingfunction ? Vars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid)) - + equation = Equation(params, grid) - + return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev) end @@ -126,7 +127,7 @@ end """ ForcedVars(dev, grid) -Returns the vars for forced two-dimensional turbulence on device dev and with +Returns the vars for forced two-dimensional turbulence on device dev and with `grid`. """ function ForcedVars(dev::Dev, grid::AbstractGrid) where Dev @@ -302,4 +303,16 @@ Returns the extraction of domain-averaged energy by drag/hypodrag μ. return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end +""" + drag_ens(prob) + +Returns the extraction of domain-averaged enstrophy by drag/hypodrag μ. +""" +@inline function drag_ens(prob) + sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid + @. vars.uh = grid.Krsq^params.nμ * abs2(sol) + vars.uh[1, 1] = 0 + return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + end # module From b583400ce5321c0f7ed94623e2573780607f9353 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Wed, 29 Jul 2020 10:29:40 -0700 Subject: [PATCH 04/19] Update twodnavierstokes.jl --- src/twodnavierstokes.jl | 35 ++++++++++++++++++++++++++++++++++- 1 file changed, 34 insertions(+), 1 deletion(-) diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl index 181f8aa0..a999f1a8 100644 --- a/src/twodnavierstokes.jl +++ b/src/twodnavierstokes.jl @@ -10,7 +10,9 @@ export dissipation, work, drag, - drag_ens + drag_ens, + work_ens, + dissipation_ens using Reexport @@ -272,6 +274,18 @@ Returns the domain-averaged dissipation rate. nν must be >= 1. return params.ν / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end +""" + dissipation_ens(prob) + +Returns the domain-averaged dissipation rate of enstrophy. nν must be >= 1. +""" +@inline function dissipation_ens(prob) + sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid + @. vars.uh = grid.Krsq^params.nν * abs2(sol) + vars.uh[1, 1] = 0 + return params.ν / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + """ work(prob) work(sol, v, grid) @@ -291,6 +305,25 @@ end @inline work(prob) = work(prob.sol, prob.vars, prob.grid) +""" + work_ens(prob) + work_ens(sol, v, grid) + +Returns the domain-averaged rate of work of enstrophy by the forcing Fh. +""" +@inline function work_ens(sol, vars::ForcedVars, grid) + @. vars.uh = sol * conj(vars.Fh) + return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + +@inline function work_ens(sol, vars::StochasticForcedVars, grid) + @. vars.uh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich + # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + +@inline work_ens(prob) = work_ens(prob.sol, prob.vars, prob.grid) + """ drag(prob) From 3325b75647c9386818eaad4d0a1809e32923db26 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Fri, 4 Sep 2020 16:43:43 -0700 Subject: [PATCH 05/19] Update barotropicqg.jl --- src/barotropicqg.jl | 58 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 52 insertions(+), 6 deletions(-) diff --git a/src/barotropicqg.jl b/src/barotropicqg.jl index 9b98e718..5d99999b 100644 --- a/src/barotropicqg.jl +++ b/src/barotropicqg.jl @@ -11,7 +11,10 @@ export meanenstrophy, dissipation, work, - drag + drag, + drag_ens, + work_ens, + dissipation_ens using Reexport @@ -60,9 +63,9 @@ function Problem(dev::Device=CPU(); params = !(typeof(eta)<:ArrayType(dev)) ? Params(grid, β, eta, μ, ν, nν, calcFU, calcFq) : Params(β, eta, rfft(eta), μ, ν, nν, calcFU, calcFq) vars = (calcFq == nothingfunction && calcFU == nothingfunction) ? Vars(dev, grid) : (stochastic ? StochasticForcedVars(dev, grid) : ForcedVars(dev, grid)) - + equation = Equation(params, grid) - + return FourierFlows.Problem(equation, stepper, dt, grid, vars, params, dev) end @@ -204,7 +207,7 @@ function calcN_advection!(N, sol, t, clock, vars, params, grid) @. vars.v = vars.v * vars.q # v*q mul!(vars.uh, grid.rfftplan, vars.u) # \hat{(u+U)*q} - + # Nonlinear advection term for q (part 1) @. N = -im * grid.kr * vars.uh # -∂[(U+u)q]/∂x mul!(vars.uh, grid.rfftplan, vars.v) # \hat{v*q} @@ -339,18 +342,30 @@ enstrophy(prob) = enstrophy(prob.sol, prob.grid, prob.vars) """ meanenergy(prob) - + Returns the energy of the domain-averaged U. """ meanenergy(prob) = real(0.5 * prob.sol[1, 1].^2) """ meanenstrophy(prob) - + Returns the enstrophy of the domain-averaged U. """ meanenstrophy(prob) = real(prob.params.β * prob.sol[1, 1]) +""" + dissipation_ens(prob) + +Returns the domain-averaged dissipation rate of enstrophy. nν must be >= 1. +""" +@inline function dissipation_ens(prob) + sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid + @. vars.uh = grid.Krsq^params.nν * abs2(sol) + vars.uh[1, 1] = 0 + return params.ν / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + """ dissipation(prob) dissipation(sol, vars, params, grid) @@ -384,6 +399,25 @@ end @inline work(prob) = work(prob.sol, prob.vars, prob.grid) +""" + work_ens(prob) + work_ens(sol, v, grid) + +Returns the domain-averaged rate of work of enstrophy by the forcing Fh. +""" +@inline function work_ens(sol, vars::ForcedVars, grid) + @. vars.uh = sol * conj(vars.Fh) + return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + +@inline function work_ens(sol, vars::StochasticForcedVars, grid) + @. vars.uh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich + # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito + return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + +@inline work_ens(prob) = work_ens(prob.sol, prob.vars, prob.grid) + """ drag(prob) drag(sol, vars, params, grid) @@ -397,4 +431,16 @@ Returns the extraction of domain-averaged energy by drag μ. return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end +""" + drag_ens(prob) + +Returns the extraction of domain-averaged enstrophy by drag/hypodrag μ. +""" +@inline function drag_ens(prob) + sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid + @. vars.uh = grid.Krsq^0.0 * abs2(sol) + vars.uh[1, 1] = 0 + return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) +end + end # module From fa44810eab8415f9e4ad9833c5d0b62bbbc261e2 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Fri, 4 Sep 2020 16:54:03 -0700 Subject: [PATCH 06/19] Update barotropicqg.jl --- src/barotropicqg.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/barotropicqg.jl b/src/barotropicqg.jl index 5d99999b..7b05c401 100644 --- a/src/barotropicqg.jl +++ b/src/barotropicqg.jl @@ -406,12 +406,12 @@ end Returns the domain-averaged rate of work of enstrophy by the forcing Fh. """ @inline function work_ens(sol, vars::ForcedVars, grid) - @. vars.uh = sol * conj(vars.Fh) + @. vars.uh = sol * conj(vars.Fqh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end @inline function work_ens(sol, vars::StochasticForcedVars, grid) - @. vars.uh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich + @. vars.uh = (vars.prevsol + sol) / 2 * conj(vars.Fqh) # Stratonovich # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end From eb2a94628bd6b8134633a970d11c83426aa4381c Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Mon, 7 Sep 2020 17:07:43 -0700 Subject: [PATCH 07/19] Tidy changes --- .gitignore | 2 +- examples/2D_Paper/2D_Paper_Inverse_Cascade.jl | 328 ------------------ 2 files changed, 1 insertion(+), 329 deletions(-) delete mode 100644 examples/2D_Paper/2D_Paper_Inverse_Cascade.jl diff --git a/.gitignore b/.gitignore index 5cd30ca1..59fc4100 100644 --- a/.gitignore +++ b/.gitignore @@ -19,6 +19,7 @@ docs/Manifest.toml # Model output **/*.jld **/*.jld2 +**/*.mat # Plots and movies **/*.png @@ -38,4 +39,3 @@ coverage/ # MacOS stuff *.DS_Store -*.mat diff --git a/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl b/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl deleted file mode 100644 index 52a77944..00000000 --- a/examples/2D_Paper/2D_Paper_Inverse_Cascade.jl +++ /dev/null @@ -1,328 +0,0 @@ -# # 2D forced-dissipative turbulence -# -#md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). -#md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). -# -# A simulation of forced-dissipative two-dimensional turbulence. We solve the -# two-dimensional vorticity equation with stochastic excitation and dissipation in -# the form of linear drag and hyperviscosity. As a demonstration, we compute how -# each of the forcing and dissipation terms contribute to the energy budget. - -using FourierFlows, Printf, Random, Plots - -using FourierFlows: parsevalsum -using Random: seed! -using FFTW: rfft, irfft - -using MAT - -import GeophysicalFlows.TwoDNavierStokes -import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy, dissipation, work, drag -import GeophysicalFlows: peakedisotropicspectrum - - -# ## Choosing a device: CPU or GPU - -dev = CPU() # Device (CPU/GPU) -nothing # hide - - -# ## Numerical, domain, and simulation parameters -# -# First, we pick some numerical and physical parameters for our model. - - n, L = 256, 2π # grid resolution and domain length - ν, nν = 2e-7, 2 # hyperviscosity coefficient and order - μ, nμ = 1e-1, 0 # linear drag coefficient -dt, tf = 0.005, 0.2/μ # timestep and final time - nt = round(Int, tf/dt) # total timesteps - ns = 4 # how many intermediate times we want to plot - nsubs = round(Int, nt/ns) -nothing # hide - - -# ## Forcing - -# We force the vorticity equation with stochastic excitation that is delta-correlated -# in time and while spatially homogeneously and isotropically correlated. The forcing -# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and -# width $\delta k_f$, and it injects energy per unit area and per unit time equal -# to $\varepsilon$. - -forcing_wavenumber = 14.0 # the central forcing wavenumber for a spectrum that is a ring in wavenumber space -forcing_bandwidth = 1.5 # the width of the forcing spectrum -ε = 0.001 # energy input rate by the forcing - -gr = TwoDGrid(dev, n, L) -x, y = gr.x, gr.y - -forcing_spectrum = @. exp(-(sqrt(gr.Krsq)-forcing_wavenumber)^2/(2*forcing_bandwidth^2)) -forcing_spectrum[ gr.Krsq .< (2π/L*2)^2 ] .= 0 # make sure that focing has no power for low wavenumbers -forcing_spectrum[ gr.Krsq .> (2π/L*20)^2 ] .= 0 # make sure that focing has no power for high wavenumbers -ε0 = parsevalsum(forcing_spectrum.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly) -forcing_spectrum .= ε/ε0 * forcing_spectrum # normalize forcing to inject energy ε - -seed!(1234) -nothing # hide - -# Next we construct function `calcF!` that computes a forcing realization every timestep -function calcF!(Fh, sol, t, cl, v, p, g) - ξ = ArrayType(dev)(exp.(2π*im*rand(eltype(gr), size(sol)))/sqrt(cl.dt)) - ξ[1, 1] = 0 - @. Fh = ξ*sqrt(forcing_spectrum) - nothing -end -nothing # hide - - -# ## Problem setup -# We initialize a `Problem` by providing a set of keyword arguments. The -# `stepper` keyword defines the time-stepper to be used. -prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, stepper="ETDRK4", - calcF=calcF!, stochastic=true) -nothing # hide - -# Define some shortcuts for convenience. -sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid -nothing # hide - - -# First let's see how a forcing realization looks like. -calcF!(v.Fh, sol, 0.0, cl, v, p, g) - -heatmap(x, y, irfft(v.Fh, g.nx), - aspectratio = 1, - c = :balance, - clim = (-200, 200), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "a forcing realization", - framestyle = :box) - - -# ## Setting initial conditions - -# Our initial condition is simply fluid at rest. -TwoDNavierStokes.set_zeta!(prob, zeros(g.nx, g.ny)) - - -# ## Diagnostics - -# Create Diagnostics; the diagnostics are aimed to probe the energy budget. -E = Diagnostic(energy, prob, nsteps=nt) # energy -R = Diagnostic(drag, prob, nsteps=nt) # dissipation by drag -D = Diagnostic(dissipation, prob, nsteps=nt) # dissipation by hyperviscosity -W = Diagnostic(work, prob, nsteps=nt) # work input by -Z = Diagnostic(enstrophy, prob; nsteps=nt) -diags = [E, D, W, R, Z] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. -nothing # hide - - -# ## Visualizing the simulation - -# We define a function that plots the vorticity field and the evolution of -# the diagnostics: energy and all terms involved in the energy budget. Last -# we confirm whether the energy budget is accurate, i.e., $\mathrm{d}E/\mathrm{d}t = W - R - D$. - -function computetendencies_and_makeplot(prob, diags) - TwoDNavierStokes.updatevars!(prob) - E, D, W, R, Z = diags - - clocktime = round(μ*cl.t, digits=2) - - i₀ = 1 - dEdt_numerical = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt #numerical first-order approximation of energy tendency - ii = (i₀):E.i-1 - ii2 = (i₀+1):E.i - - t = E.t[ii] - dEdt_computed = W[ii2] - D[ii] - R[ii] # Stratonovich interpretation - - residual = dEdt_computed - dEdt_numerical - - l = @layout grid(2, 2) - - p1 = heatmap(x, y, v.zeta, - aspectratio = 1, - legend = false, - c = :viridis, - clim = (-25, 25), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "μt", - ylabel = "y", - title = "∇²ψ(x, y, t="*@sprintf("%.2f", cl.t)*")", - framestyle = :box) - - p2 = plot(μ*t, [W[ii2] ε.+0*t -D[ii] -R[ii]], - label = ["work, W" "ensemble mean work, " "dissipation, D" "drag, D=-2μE"], - linestyle = [:solid :dash :solid :solid], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "energy sources and sinks") - - p3 = plot(μ*t, [dEdt_computed[ii], dEdt_numerical], - label = ["computed W-D" "numerical dE/dt"], - linestyle = [:solid :dashdotdot], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "dE/dt") - - p4 = plot(μ*t, residual, - label = "residual dE/dt = computed - numerical", - linewidth = 2, - alpha = 0.7, - xlabel = "μt") - - p = plot(p1, p2, p3, p4, layout=l, size = (900, 800)) - return p -end -nothing # hide - - - - -# ## Time-stepping the `Problem` forward - -# Finally, we time-step the `Problem` forward in time. - -startwalltime = time() -for i = 1:ns - stepforward!(prob, diags, round(Int, nt/ns)) - TwoDNavierStokes.updatevars!(prob) - cfl = cl.dt*maximum([maximum(v.u)/g.dx, maximum(v.v)/g.dy]) - - log = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min", cl.step, cl.t, - cfl, (time()-startwalltime)/60) - - println(log) -end - - -# ## Plot -# And now let's see what we got. We plot the output. - -p = computetendencies_and_makeplot(prob, diags) - -# ## Output - -# We choose folder for outputing `.jld2` files and snapshots (`.png` files). -filepath = "." -plotpath = "/Users/brodiepearson/GitHub/GeophysicalFlows.jl/examples/2D_Paper/Plots/Inverse_Cascade" -plotname = "snapshots" -filename = joinpath(filepath, "Inverse_Cascade.jld2") -matfilename = joinpath(filepath, "Inverse_Cascade.mat") -nothing # hide - -# Do some basic file management -if isfile(filename); rm(filename); end -if !isdir(plotpath); mkdir(plotpath); end -nothing # hide - -# And then create Output -get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution -get_u(prob) = Array(irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx)) -out = Output(prob, filename, (:sol, get_sol), (:u, get_u)) -saveproblem(out) -matfile = matopen(matfilename) -u_temp = (:u, get_u) -write(matfile, "zeta", v.zeta) -write(matfile, "u", v.u) -write(matfile, "v", v.v) -write(matfile, "nx", gr.nx) -write(matfile, "ny", gr.ny) -write(matfile, "dx", gr.dx) -write(matfile, "dy", gr.dy) -write(matfile, "Lx", gr.Lx) -write(matfile, "Ly", gr.Ly) -close(matfile) -nothing # hide - -# ## Visualizing the simulation - -# We initialize a plot with the vorticity field and the time-series of -# energy and enstrophy diagnostics. - -p1 = heatmap(x, y, v.zeta, - aspectratio = 1, - c = :balance, - clim = (-40, 40), - xlims = (-L/2, L/2), - ylims = (-L/2, L/2), - xticks = -3:3, - yticks = -3:3, - xlabel = "x", - ylabel = "y", - title = "vorticity, t="*@sprintf("%.2f", cl.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 grid(1, 2) -p = plot(p1, p2, layout = l, size = (900, 400)) - - -# ## Time-stepping the `Problem` forward - -# We time-step the `Problem` forward in time. - -startwalltime = time() - -anim = @animate for j = 0:Int(nt/nsubs) - - log = @sprintf("step: %04d, t: %d, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min", - cl.step, cl.t, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60) - - if j%(1000/nsubs)==0; println(log) end - - p[1][1][:z] = v.zeta - p[1][:title] = "vorticity, t="*@sprintf("%.2f", cl.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]) - - stepforward!(prob, diags, nsubs) - TwoDNavierStokes.updatevars!(prob) - -end - -mp4(anim, "twodturb.mp4", fps=18) - - -# Last we save the output. -saveoutput(out) - - -# ## Radial energy spectrum - -# After the simulation is done we plot the radial energy spectrum to illustrate -# how `FourierFlows.radialspectrum` can be used, - -E = @. 0.5*(v.u^2 + v.v^2) # energy density -Eh = rfft(E) # Fourier transform of energy density -kr, Ehr = FourierFlows.radialspectrum(Eh, gr, refinement=1) # compute radial specturm of `Eh` -nothing # hide - -# and we plot it. -plot(kr, abs.(Ehr), - linewidth = 2, - alpha = 0.7, - xlabel = "kᵣ", ylabel = "∫ |Ê| kᵣ dk_θ", - xlims = (5e-1, gr.nx), - xscale = :log10, yscale = :log10, - title = "Radial energy spectrum", - legend = false) From c4c9436f49d1c046a9ff78329b5c080368dca4d6 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Tue, 8 Sep 2020 23:44:04 -0700 Subject: [PATCH 08/19] Add Enstrophy Budget Analyses to Forced 2D Example --- .../twodnavierstokes_stochasticforcing.jl | 99 +++++++++++++------ 1 file changed, 68 insertions(+), 31 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 2541e4c5..5bfa5f7a 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -3,19 +3,20 @@ #md # This example can be run online via [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). #md # Also, it can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/twodnavierstokes_stochasticforcing.ipynb). # -# A simulation of forced-dissipative two-dimensional turbulence. We solve the +# A simulation of forced-dissipative two-dimensional turbulence. We solve the # two-dimensional vorticity equation with stochastic excitation and dissipation in -# the form of linear drag and hyperviscosity. As a demonstration, we compute how +# the form of linear drag and hyperviscosity. As a demonstration, we compute how # each of the forcing and dissipation terms contribute to the energy budget. using FourierFlows, Printf, Plots - + using FourierFlows: parsevalsum using Random: seed! using FFTW: irfft import GeophysicalFlows.TwoDNavierStokes import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy, dissipation, work, drag +import GeophysicalFlows.TwoDNavierStokes: dissipation_ens, work_ens, drag_ens # ## Choosing a device: CPU or GPU @@ -41,12 +42,12 @@ nothing # hide # We force the vorticity equation with stochastic excitation that is delta-correlated # in time and while spatially homogeneously and isotropically correlated. The forcing -# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and -# width $\delta k_f$, and it injects energy per unit area and per unit time equal +# has a spectrum with power in a ring in wavenumber space of radious $k_f$ and +# width $\delta k_f$, and it injects energy per unit area and per unit time equal # to $\varepsilon$. forcing_wavenumber = 14.0 # the central forcing wavenumber for a spectrum that is a ring in wavenumber space -forcing_bandwidth = 1.5 # the width of the forcing spectrum +forcing_bandwidth = 1.5 # the width of the forcing spectrum ε = 0.001 # energy input rate by the forcing gr = TwoDGrid(dev, n, L) @@ -113,7 +114,11 @@ E = Diagnostic(energy, prob, nsteps=nt) # energy R = Diagnostic(drag, prob, nsteps=nt) # dissipation by drag D = Diagnostic(dissipation, prob, nsteps=nt) # dissipation by hyperviscosity W = Diagnostic(work, prob, nsteps=nt) # work input by forcing -diags = [E, D, W, R] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. +Z = Diagnostic(enstrophy, prob, nsteps=nt) # energy +RZ = Diagnostic(drag_ens, prob, nsteps=nt) # dissipation by drag +DZ = Diagnostic(dissipation_ens, prob, nsteps=nt) # dissipation by hyperviscosity +WZ = Diagnostic(work_ens, prob, nsteps=nt) # work input by forcing +diags = [E, D, W, R, Z, DZ, WZ, RZ] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -125,23 +130,28 @@ nothing # hide function computetendencies_and_makeplot(prob, diags) TwoDNavierStokes.updatevars!(prob) - E, D, W, R = diags + E, D, W, R, Z, DZ, WZ, RZ = diags clocktime = round(μ*cl.t, digits=2) i₀ = 1 dEdt_numerical = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt #numerical first-order approximation of energy tendency + dZdt_numerical = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt #numerical first-order approximation of enstrophy tendency ii = (i₀):E.i-1 ii2 = (i₀+1):E.i t = E.t[ii] dEdt_computed = W[ii2] - D[ii] - R[ii] # Stratonovich interpretation - - residual = dEdt_computed - dEdt_numerical + dZdt_computed = WZ[ii2] - DZ[ii] - RZ[ii] + + residual_E = dEdt_computed - dEdt_numerical + residual_Z = dZdt_computed - dZdt_numerical + + εZ = ε*(forcing_wavenumber^2) # Estimate of mean enstrophy injection rate - l = @layout grid(2, 2) + l = @layout grid(2,3) - p1 = heatmap(x, y, v.zeta, + pzeta = heatmap(x, y, v.zeta, aspectratio = 1, legend = false, c = :viridis, @@ -155,7 +165,9 @@ function computetendencies_and_makeplot(prob, diags) title = "∇²ψ(x, y, t="*@sprintf("%.2f", cl.t)*")", framestyle = :box) - p2 = plot(μ*t, [W[ii2] ε.+0*t -D[ii] -R[ii]], + pζ = plot(pzeta, size = (400, 400)) + + p1 = plot(μ*t, [W[ii2] ε.+0*t -D[ii] -R[ii]], label = ["work, W" "ensemble mean work, " "dissipation, D" "drag, D=-2μE"], linestyle = [:solid :dash :solid :solid], linewidth = 2, @@ -163,22 +175,47 @@ function computetendencies_and_makeplot(prob, diags) xlabel = "μt", ylabel = "energy sources and sinks") - p3 = plot(μ*t, [dEdt_computed[ii], dEdt_numerical], - label = ["computed W-D" "numerical dE/dt"], - linestyle = [:solid :dashdotdot], - linewidth = 2, - alpha = 0.8, - xlabel = "μt", - ylabel = "dE/dt") - - p4 = plot(μ*t, residual, - label = "residual dE/dt = computed - numerical", - linewidth = 2, - alpha = 0.7, - xlabel = "μt") - - p = plot(p1, p2, p3, p4, layout=l, size = (900, 800)) - return p + p2 = plot(μ*t, [dEdt_computed[ii], dEdt_numerical], + label = ["computed W-D" "numerical dE/dt"], + linestyle = [:solid :dashdotdot], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "dE/dt") + + p3 = plot(μ*t, residual_E, + label = "residual dE/dt = computed - numerical", + linewidth = 2, + alpha = 0.7, + xlabel = "μt") + + p4 = plot(μ*t, [WZ[ii2] εZ.+0*t -DZ[ii] -RZ[ii]], + label = ["Enstrophy work, WZ" "mean enstrophy work, " "enstrophy dissipation, DZ" "enstrophy drag, D=-2μZ"], + linestyle = [:solid :dash :solid :solid], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "enstrophy sources and sinks") + + + p5 = plot(μ*t, [dZdt_computed[ii], dZdt_numerical], + label = ["computed WZ-DZ" "numerical dZ/dt"], + linestyle = [:solid :dashdotdot], + linewidth = 2, + alpha = 0.8, + xlabel = "μt", + ylabel = "dZ/dt") + + p6 = plot(μ*t, residual_Z, + label = "residual dZ/dt = computed - numerical", + linewidth = 2, + alpha = 0.7, + xlabel = "μt") + + + pbudgets = plot(p1, p2, p3, p4, p5, p6, layout=l, size = (1300, 900)) + + return pζ, pbudgets end nothing # hide @@ -198,11 +235,11 @@ for i = 1:ns log = @sprintf("step: %04d, t: %.1f, cfl: %.3f, walltime: %.2f min", cl.step, cl.t, cfl, (time()-startwalltime)/60) - println(log) + println(log) end # ## Plot # And now let's see what we got. We plot the output. -p = computetendencies_and_makeplot(prob, diags) +pζ, pbudgets = computetendencies_and_makeplot(prob, diags) From a355eceee6861c41cc3199e331937e02f2f9781b Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Tue, 8 Sep 2020 23:44:51 -0700 Subject: [PATCH 09/19] Update Project.toml --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index ff2db8c2..06639937 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" From cb11ee424f2363d66d44b5f160db95e1243a6343 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Wed, 9 Sep 2020 00:10:00 -0700 Subject: [PATCH 10/19] Change diagnostic names in 2D model --- src/twodnavierstokes.jl | 48 ++++++++++++++++++++--------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl index a999f1a8..658244f2 100644 --- a/src/twodnavierstokes.jl +++ b/src/twodnavierstokes.jl @@ -6,13 +6,13 @@ export updatevars!, energy, + energy_dissipation, + energy_work, + energy_drag, enstrophy, - dissipation, - work, - drag, - drag_ens, - work_ens, - dissipation_ens + enstrophy_dissipation, + enstrophy_work, + enstrophy_drag using Reexport @@ -263,11 +263,11 @@ solution `sol`. end """ - dissipation(prob) + energy_dissipation(prob) Returns the domain-averaged dissipation rate. nν must be >= 1. """ -@inline function dissipation(prob) +@inline function energy_dissipation(prob) sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid @. vars.uh = grid.Krsq^(params.nν-1) * abs2(sol) vars.uh[1, 1] = 0 @@ -275,11 +275,11 @@ Returns the domain-averaged dissipation rate. nν must be >= 1. end """ - dissipation_ens(prob) + enstrophy_dissipation(prob) Returns the domain-averaged dissipation rate of enstrophy. nν must be >= 1. """ -@inline function dissipation_ens(prob) +@inline function enstrophy_dissipation(prob) sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid @. vars.uh = grid.Krsq^params.nν * abs2(sol) vars.uh[1, 1] = 0 @@ -287,49 +287,49 @@ Returns the domain-averaged dissipation rate of enstrophy. nν must be >= 1. end """ - work(prob) - work(sol, v, grid) + energy_work(prob) + energy_work(sol, v, grid) Returns the domain-averaged rate of work of energy by the forcing Fh. """ -@inline function work(sol, vars::ForcedVars, grid) +@inline function energy_work(sol, vars::ForcedVars, grid) @. vars.uh = grid.invKrsq * sol * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end -@inline function work(sol, vars::StochasticForcedVars, grid) +@inline function energy_work(sol, vars::StochasticForcedVars, grid) @. vars.uh = grid.invKrsq * (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end -@inline work(prob) = work(prob.sol, prob.vars, prob.grid) +@inline energy_work(prob) = energy_work(prob.sol, prob.vars, prob.grid) """ - work_ens(prob) - work_ens(sol, v, grid) + enstrophy_work(prob) + enstrophy_work(sol, v, grid) Returns the domain-averaged rate of work of enstrophy by the forcing Fh. """ -@inline function work_ens(sol, vars::ForcedVars, grid) +@inline function enstrophy_work(sol, vars::ForcedVars, grid) @. vars.uh = sol * conj(vars.Fh) return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end -@inline function work_ens(sol, vars::StochasticForcedVars, grid) +@inline function enstrophy_work(sol, vars::StochasticForcedVars, grid) @. vars.uh = (vars.prevsol + sol) / 2 * conj(vars.Fh) # Stratonovich # @. vars.uh = grid.invKrsq * vars.prevsol * conj(vars.Fh) # Ito return 1 / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end -@inline work_ens(prob) = work_ens(prob.sol, prob.vars, prob.grid) +@inline enstrophy_work(prob) = enstrophy_work(prob.sol, prob.vars, prob.grid) """ - drag(prob) + energy_drag(prob) Returns the extraction of domain-averaged energy by drag/hypodrag μ. """ -@inline function drag(prob) +@inline function energy_drag(prob) sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid @. vars.uh = grid.Krsq^(params.nμ-1) * abs2(sol) vars.uh[1, 1] = 0 @@ -337,11 +337,11 @@ Returns the extraction of domain-averaged energy by drag/hypodrag μ. end """ - drag_ens(prob) + enstrophy_drag(prob) Returns the extraction of domain-averaged enstrophy by drag/hypodrag μ. """ -@inline function drag_ens(prob) +@inline function enstrophy_drag(prob) sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid @. vars.uh = grid.Krsq^params.nμ * abs2(sol) vars.uh[1, 1] = 0 From 23dbd0fe701f895bd59c1c3ccb98773a97bcf235 Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Wed, 9 Sep 2020 00:13:41 -0700 Subject: [PATCH 11/19] Add enstrophy budget diagnostics --- examples/twodnavierstokes_stochasticforcing.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 5bfa5f7a..9aa78c0e 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -15,8 +15,8 @@ using Random: seed! using FFTW: irfft import GeophysicalFlows.TwoDNavierStokes -import GeophysicalFlows.TwoDNavierStokes: energy, enstrophy, dissipation, work, drag -import GeophysicalFlows.TwoDNavierStokes: dissipation_ens, work_ens, drag_ens +import GeophysicalFlows.TwoDNavierStokes: energy, energy_dissipation, energy_work, energy_drag +import GeophysicalFlows.TwoDNavierStokes: enstrophy, enstrophy_dissipation, enstrophy_work, enstrophy_drag # ## Choosing a device: CPU or GPU @@ -111,13 +111,13 @@ TwoDNavierStokes.set_zeta!(prob, zeros(g.nx, g.ny)) # Create Diagnostics; the diagnostics are aimed to probe the energy budget. E = Diagnostic(energy, prob, nsteps=nt) # energy -R = Diagnostic(drag, prob, nsteps=nt) # dissipation by drag -D = Diagnostic(dissipation, prob, nsteps=nt) # dissipation by hyperviscosity -W = Diagnostic(work, prob, nsteps=nt) # work input by forcing +R = Diagnostic(energy_drag, prob, nsteps=nt) # dissipation by drag +D = Diagnostic(energy_dissipation, prob, nsteps=nt) # dissipation by hyperviscosity +W = Diagnostic(energy_work, prob, nsteps=nt) # work input by forcing Z = Diagnostic(enstrophy, prob, nsteps=nt) # energy -RZ = Diagnostic(drag_ens, prob, nsteps=nt) # dissipation by drag -DZ = Diagnostic(dissipation_ens, prob, nsteps=nt) # dissipation by hyperviscosity -WZ = Diagnostic(work_ens, prob, nsteps=nt) # work input by forcing +RZ = Diagnostic(enstrophy_drag, prob, nsteps=nt) # dissipation by drag +DZ = Diagnostic(enstrophy_dissipation, prob, nsteps=nt) # dissipation by hyperviscosity +WZ = Diagnostic(enstrophy_work, prob, nsteps=nt) # work input by forcing diags = [E, D, W, R, Z, DZ, WZ, RZ] # A list of Diagnostics types passed to "stepforward!" will be updated every timestep. nothing # hide @@ -243,3 +243,5 @@ end # And now let's see what we got. We plot the output. pζ, pbudgets = computetendencies_and_makeplot(prob, diags) + +pbudgets From e3bb62e47cb3e9be1b17f43e279f345b0ed54f6c Mon Sep 17 00:00:00 2001 From: Brodie Pearson Date: Wed, 9 Sep 2020 00:14:14 -0700 Subject: [PATCH 12/19] Add tests for residual in enstrophy budgets --- test/test_twodnavierstokes.jl | 43 ++++++++++++++++++++++++----------- 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/test/test_twodnavierstokes.jl b/test/test_twodnavierstokes.jl index 35e5c7ff..611051dd 100644 --- a/test/test_twodnavierstokes.jl +++ b/test/test_twodnavierstokes.jl @@ -56,19 +56,24 @@ function test_twodnavierstokes_stochasticforcingbudgets(dev::Device=CPU(); n=256 TwoDNavierStokes.set_zeta!(prob, 0*x) E = Diagnostic(TwoDNavierStokes.energy, prob, nsteps=nt) - D = Diagnostic(TwoDNavierStokes.dissipation, prob, nsteps=nt) - R = Diagnostic(TwoDNavierStokes.drag, prob, nsteps=nt) - W = Diagnostic(TwoDNavierStokes.work, prob, nsteps=nt) - diags = [E, D, W, R] + D = Diagnostic(TwoDNavierStokes.energy_dissipation, prob, nsteps=nt) + R = Diagnostic(TwoDNavierStokes.energy_drag, prob, nsteps=nt) + W = Diagnostic(TwoDNavierStokes.energy_work, prob, nsteps=nt) + Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nt) + DZ = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt) + RZ = Diagnostic(TwoDNavierStokes.enstrophy_drag, prob, nsteps=nt) + WZ = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, nsteps=nt) + diags = [E, D, W, R, Z, DZ, WZ, RZ] stepforward!(prob, diags, nt) TwoDNavierStokes.updatevars!(prob) - E, D, W, R = diags + E, D, W, R, Z, DZ, WZ, RZ = diags t = round(μ*cl.t, digits=2) i₀ = 1 dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt + dZdt = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt ii = (i₀):E.i-1 ii2 = (i₀+1):E.i @@ -77,9 +82,13 @@ function test_twodnavierstokes_stochasticforcingbudgets(dev::Device=CPU(); n=256 # then we need to add the drift term # total = W[ii2]+ε - D[ii] - R[ii] # Ito total = W[ii2] - D[ii] - R[ii] # Stratonovich + totalZ = WZ[ii2] - DZ[ii] - RZ[ii] residual = dEdt - total isapprox(mean(abs.(residual)), 0, atol=1e-4) + + residualZ = dZdt - totalZ + isapprox(mean(abs.(residualZ)), 0, atol=(kf^2)*1e-4) end @@ -89,7 +98,7 @@ function test_twodnavierstokes_deterministicforcingbudgets(dev::Device=CPU(); n= μ, nμ = 1e-1, 0 dt, tf = 0.005, 0.1/μ nt = round(Int, tf/dt) - + gr = TwoDGrid(dev, n, L) x, y = gridpoints(gr) @@ -109,28 +118,36 @@ function test_twodnavierstokes_deterministicforcingbudgets(dev::Device=CPU(); n= TwoDNavierStokes.set_zeta!(prob, 0*x) E = Diagnostic(TwoDNavierStokes.energy, prob, nsteps=nt) - D = Diagnostic(TwoDNavierStokes.dissipation, prob, nsteps=nt) - R = Diagnostic(TwoDNavierStokes.drag, prob, nsteps=nt) - W = Diagnostic(TwoDNavierStokes.work, prob, nsteps=nt) - diags = [E, D, W, R] + D = Diagnostic(TwoDNavierStokes.energy_dissipation, prob, nsteps=nt) + R = Diagnostic(TwoDNavierStokes.energy_drag, prob, nsteps=nt) + W = Diagnostic(TwoDNavierStokes.energy_work, prob, nsteps=nt) + Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nt) + DZ = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt) + RZ = Diagnostic(TwoDNavierStokes.enstrophy_drag, prob, nsteps=nt) + WZ = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, nsteps=nt) + diags = [E, D, W, R, Z, DZ, WZ, RZ] # Step forward stepforward!(prob, diags, nt) TwoDNavierStokes.updatevars!(prob) - E, D, W, R = diags + E, D, W, R, Z, DZ, WZ, RZ = diags t = round(μ*cl.t, digits=2) i₀ = 1 dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt + dZdt = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt ii = (i₀):E.i-1 ii2 = (i₀+1):E.i # dEdt = W - D - R? total = W[ii2] - D[ii] - R[ii] + totalZ = WZ[ii2] - DZ[ii] - RZ[ii] residual = dEdt - total + residualZ = dZdt - totalZ isapprox(mean(abs.(residual)), 0, atol=1e-8) + isapprox(mean(abs.(residualZ)), 0, atol=1e-8) end """ @@ -194,7 +211,7 @@ function test_twodnavierstokes_energyenstrophy(dev::Device=CPU()) enstrophy_calc = 2701/162 prob = TwoDNavierStokes.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, stepper="ForwardEuler") - + sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid; TwoDNavierStokes.set_zeta!(prob, ζ₀) @@ -202,7 +219,7 @@ function test_twodnavierstokes_energyenstrophy(dev::Device=CPU()) energyζ₀ = TwoDNavierStokes.energy(prob) enstrophyζ₀ = TwoDNavierStokes.enstrophy(prob) - + params = TwoDNavierStokes.Params(p.ν, p.nν) (isapprox(energyζ₀, energy_calc, rtol=rtol_twodnavierstokes) && From 898ffb158d84cd06a7c30be420c21aceb0f145e5 Mon Sep 17 00:00:00 2001 From: BrodieP Date: Thu, 10 Sep 2020 22:27:22 -0700 Subject: [PATCH 13/19] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index dfb19a16..03f39ce8 100644 --- a/Project.toml +++ b/Project.toml @@ -10,7 +10,6 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" FourierFlows = "2aec4490-903f-5c70-9b11-9bed06a700e1" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" From 65449024e9b2c62b8a930f100ddc143943de363a Mon Sep 17 00:00:00 2001 From: BrodieP Date: Thu, 10 Sep 2020 22:39:18 -0700 Subject: [PATCH 14/19] Update twodnavierstokes_stochasticforcing.jl Added more accurate enstrophy dissipation. Fixed legends to fix typos and remove ambiguity. Changed enstrophy diagnostics to superscript Z notation (otherwise the drag legend was ambiguous) --- examples/twodnavierstokes_stochasticforcing.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index 9aa78c0e..bcb8c455 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -130,7 +130,7 @@ nothing # hide function computetendencies_and_makeplot(prob, diags) TwoDNavierStokes.updatevars!(prob) - E, D, W, R, Z, DZ, WZ, RZ = diags + E, D, W, R, Z, Dᶻ, Wᶻ, Rᶻ = diags clocktime = round(μ*cl.t, digits=2) @@ -142,12 +142,12 @@ function computetendencies_and_makeplot(prob, diags) t = E.t[ii] dEdt_computed = W[ii2] - D[ii] - R[ii] # Stratonovich interpretation - dZdt_computed = WZ[ii2] - DZ[ii] - RZ[ii] + dZdt_computed = Wᶻ[ii2] - Dᶻ[ii] - Rᶻ[ii] residual_E = dEdt_computed - dEdt_numerical residual_Z = dZdt_computed - dZdt_numerical - εZ = ε*(forcing_wavenumber^2) # Estimate of mean enstrophy injection rate + εᶻ = parsevalsum(forcing_spectrum / 2, g) / (g.Lx * g.Ly) l = @layout grid(2,3) @@ -168,7 +168,7 @@ function computetendencies_and_makeplot(prob, diags) pζ = plot(pzeta, size = (400, 400)) p1 = plot(μ*t, [W[ii2] ε.+0*t -D[ii] -R[ii]], - label = ["work, W" "ensemble mean work, " "dissipation, D" "drag, D=-2μE"], + label = ["work, W" "ensemble mean work, " "dissipation, D" "drag, R=-2μE"], linestyle = [:solid :dash :solid :solid], linewidth = 2, alpha = 0.8, @@ -189,8 +189,8 @@ function computetendencies_and_makeplot(prob, diags) alpha = 0.7, xlabel = "μt") - p4 = plot(μ*t, [WZ[ii2] εZ.+0*t -DZ[ii] -RZ[ii]], - label = ["Enstrophy work, WZ" "mean enstrophy work, " "enstrophy dissipation, DZ" "enstrophy drag, D=-2μZ"], + p4 = plot(μ*t, [Wᶻ[ii2] εᶻ.+0*t -Dᶻ[ii] -Rᶻ[ii]], + label = ["Enstrophy work, Wᶻ" "mean enstrophy work, " "enstrophy dissipation, Dᶻ" "enstrophy drag, Rᶻ=-2μZ"], linestyle = [:solid :dash :solid :solid], linewidth = 2, alpha = 0.8, @@ -199,7 +199,7 @@ function computetendencies_and_makeplot(prob, diags) p5 = plot(μ*t, [dZdt_computed[ii], dZdt_numerical], - label = ["computed WZ-DZ" "numerical dZ/dt"], + label = ["computed Wᶻ-Dᶻ" "numerical dZ/dt"], linestyle = [:solid :dashdotdot], linewidth = 2, alpha = 0.8, From 0d1b8cbc74de65bf239e11e44574b3969f790c6a Mon Sep 17 00:00:00 2001 From: BrodieP Date: Thu, 10 Sep 2020 22:42:45 -0700 Subject: [PATCH 15/19] Update barotropicqg.jl Added missing CUDA statement --- src/barotropicqg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/barotropicqg.jl b/src/barotropicqg.jl index e5848903..4bde49d2 100644 --- a/src/barotropicqg.jl +++ b/src/barotropicqg.jl @@ -441,7 +441,7 @@ Returns the extraction of domain-averaged enstrophy by drag/hypodrag μ. @inline function drag_ens(prob) sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid @. vars.uh = grid.Krsq^0.0 * abs2(sol) - vars.uh[1, 1] = 0 + CUDA.@allowscalar vars.uh[1, 1] = 0 return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end From df34179e48fb7e755d0ae226dddbd629c077184c Mon Sep 17 00:00:00 2001 From: BrodieP Date: Thu, 10 Sep 2020 22:43:35 -0700 Subject: [PATCH 16/19] Update twodnavierstokes.jl Added missing CUDA statement --- src/twodnavierstokes.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/twodnavierstokes.jl b/src/twodnavierstokes.jl index 807a5175..1fe2b504 100644 --- a/src/twodnavierstokes.jl +++ b/src/twodnavierstokes.jl @@ -346,7 +346,7 @@ Returns the extraction of domain-averaged enstrophy by drag/hypodrag μ. @inline function enstrophy_drag(prob) sol, vars, params, grid = prob.sol, prob.vars, prob.params, prob.grid @. vars.uh = grid.Krsq^params.nμ * abs2(sol) - vars.uh[1, 1] = 0 + CUDA.@allowscalar vars.uh[1, 1] = 0 return params.μ / (grid.Lx * grid.Ly) * parsevalsum(vars.uh, grid) end From 5e666e1faba3a874d777319b804a4b9cdbfbc659 Mon Sep 17 00:00:00 2001 From: BrodieP Date: Thu, 10 Sep 2020 22:56:58 -0700 Subject: [PATCH 17/19] Update test_twodnavierstokes.jl Modification so check tests both energy and enstrophy budgets --- test/test_twodnavierstokes.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/test/test_twodnavierstokes.jl b/test/test_twodnavierstokes.jl index 8849e1c8..9fbdd327 100644 --- a/test/test_twodnavierstokes.jl +++ b/test/test_twodnavierstokes.jl @@ -85,10 +85,9 @@ function test_twodnavierstokes_stochasticforcingbudgets(dev::Device=CPU(); n=256 totalZ = WZ[ii2] - DZ[ii] - RZ[ii] residual = dEdt - total - isapprox(mean(abs.(residual)), 0, atol=1e-4) - residualZ = dZdt - totalZ - isapprox(mean(abs.(residualZ)), 0, atol=(kf^2)*1e-4) + + return isapprox(mean(abs.(residual)), 0, atol=1e-4) && isapprox(mean(abs.(residualZ)), 0, atol=(kf^2)*1e-4) end @@ -146,8 +145,7 @@ function test_twodnavierstokes_deterministicforcingbudgets(dev::Device=CPU(); n= residual = dEdt - total residualZ = dZdt - totalZ - isapprox(mean(abs.(residual)), 0, atol=1e-8) - isapprox(mean(abs.(residualZ)), 0, atol=1e-8) + return isapprox(mean(abs.(residual)), 0, atol=1e-8) && isapprox(mean(abs.(residualZ)), 0, atol=1e-8) end """ @@ -231,6 +229,6 @@ function test_twodnavierstokes_problemtype(dev, T) prob = TwoDNavierStokes.Problem(dev; T=T) A = ArrayType(dev) - + (typeof(prob.sol)<:A{Complex{T},2} && typeof(prob.grid.Lx)==T && eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T,2}) end From 8b357cc486a33e58cb50d1ff4a754d4631a947ef Mon Sep 17 00:00:00 2001 From: Navid Constantinou Date: Fri, 11 Sep 2020 20:29:37 +1000 Subject: [PATCH 18/19] splits enenergy and enstrophy budget into different tests --- test/runtests.jl | 6 +- test/test_twodnavierstokes.jl | 142 ++++++++++++++++++++++++++++------ 2 files changed, 122 insertions(+), 26 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 14962921..a818e11b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -48,8 +48,10 @@ for dev in devices @test test_twodnavierstokes_advection(0.0005, "ForwardEuler", dev) @test test_twodnavierstokes_lambdipole(256, 1e-3, dev) - @test test_twodnavierstokes_stochasticforcingbudgets(dev) - @test test_twodnavierstokes_deterministicforcingbudgets(dev) + @test test_twodnavierstokes_deterministicforcing_energybudget(dev) + @test test_twodnavierstokes_stochasticforcing_energybudget(dev) + @test test_twodnavierstokes_deterministicforcing_enstrophybudget(dev) + @test test_twodnavierstokes_stochasticforcing_enstrophybudget(dev) @test test_twodnavierstokes_energyenstrophy(dev) @test test_twodnavierstokes_problemtype(dev, Float32) @test TwoDNavierStokes.nothingfunction() == nothing diff --git a/test/test_twodnavierstokes.jl b/test/test_twodnavierstokes.jl index 9fbdd327..52aaa4de 100644 --- a/test/test_twodnavierstokes.jl +++ b/test/test_twodnavierstokes.jl @@ -20,7 +20,7 @@ function test_twodnavierstokes_lambdipole(n, dt, dev::Device=CPU(); L=2π, Ue=1, isapprox(Ue, mean(Ue_m[2:end]), rtol=rtol_lambdipole) end -function test_twodnavierstokes_stochasticforcingbudgets(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, μ=1e-1, nμ=0, tf=0.1/μ) +function test_twodnavierstokes_stochasticforcing_energybudget(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, μ=1e-1, nμ=0, tf=0.1/μ) nt = round(Int, tf/dt) # Forcing parameters @@ -59,21 +59,15 @@ function test_twodnavierstokes_stochasticforcingbudgets(dev::Device=CPU(); n=256 D = Diagnostic(TwoDNavierStokes.energy_dissipation, prob, nsteps=nt) R = Diagnostic(TwoDNavierStokes.energy_drag, prob, nsteps=nt) W = Diagnostic(TwoDNavierStokes.energy_work, prob, nsteps=nt) - Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nt) - DZ = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt) - RZ = Diagnostic(TwoDNavierStokes.enstrophy_drag, prob, nsteps=nt) - WZ = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, nsteps=nt) - diags = [E, D, W, R, Z, DZ, WZ, RZ] + diags = [E, D, W, R] stepforward!(prob, diags, nt) TwoDNavierStokes.updatevars!(prob) - E, D, W, R, Z, DZ, WZ, RZ = diags - t = round(μ*cl.t, digits=2) + E, D, W, R = diags i₀ = 1 dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt - dZdt = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt ii = (i₀):E.i-1 ii2 = (i₀+1):E.i @@ -82,16 +76,75 @@ function test_twodnavierstokes_stochasticforcingbudgets(dev::Device=CPU(); n=256 # then we need to add the drift term # total = W[ii2]+ε - D[ii] - R[ii] # Ito total = W[ii2] - D[ii] - R[ii] # Stratonovich - totalZ = WZ[ii2] - DZ[ii] - RZ[ii] residual = dEdt - total - residualZ = dZdt - totalZ - return isapprox(mean(abs.(residual)), 0, atol=1e-4) && isapprox(mean(abs.(residualZ)), 0, atol=(kf^2)*1e-4) + return isapprox(mean(abs.(residual)), 0, atol=1e-4) end +function test_twodnavierstokes_stochasticforcing_enstrophybudget(dev::Device=CPU(); n=256, L=2π, dt=0.005, ν=1e-7, nν=2, μ=1e-1, nμ=0, tf=0.1/μ) + nt = round(Int, tf/dt) + + # Forcing parameters + kf, dkf = 12.0, 2.0 + ε = 0.1 + + gr = TwoDGrid(dev, n, L) + x, y = gridpoints(gr) + + Kr = ArrayType(dev)([CUDA.@allowscalar gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]) + + forcingcovariancespectrum = ArrayType(dev)(zero(gr.Krsq)) + @. forcingcovariancespectrum = exp(-(sqrt(gr.Krsq) - kf)^2 / (2 * dkf^2)) + CUDA.@allowscalar @. forcingcovariancespectrum[gr.Krsq .< 2^2] = 0 + CUDA.@allowscalar @. forcingcovariancespectrum[gr.Krsq .> 20^2] = 0 + CUDA.@allowscalar @. forcingcovariancespectrum[Kr .< 2π/L] = 0 + ε0 = parsevalsum(forcingcovariancespectrum .* gr.invKrsq / 2, gr) / (gr.Lx * gr.Ly) + forcingcovariancespectrum .= ε / ε0 * forcingcovariancespectrum + + Random.seed!(1234) + + function calcF!(Fh, sol, t, cl, v, p, g) + eta = ArrayType(dev)(exp.(2π * im * rand(Float64, size(sol))) / sqrt(cl.dt)) + CUDA.@allowscalar eta[1, 1] = 0.0 + @. Fh = eta * sqrt(forcingcovariancespectrum) + nothing + end + + prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, + stepper="RK4", calcF=calcF!, stochastic=true) + + sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid; + + TwoDNavierStokes.set_zeta!(prob, 0*x) + Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nt) + D = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt) + R = Diagnostic(TwoDNavierStokes.enstrophy_drag, prob, nsteps=nt) + W = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, nsteps=nt) + diags = [Z, D, W, R] + + stepforward!(prob, diags, nt) + TwoDNavierStokes.updatevars!(prob) + + Z, D, W, R = diags + + i₀ = 1 + dZdt = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt + ii = (i₀):Z.i-1 + ii2 = (i₀+1):Z.i + + # dZdt = W - D - R? + # If the Ito interpretation was used for the work + # then we need to add the drift term + # total = W[ii2]+ε - D[ii] - R[ii] # Ito + total = W[ii2] - D[ii] - R[ii] # Stratonovich + + residual = dZdt - total + + return isapprox(mean(abs.(residual)), 0, atol=(kf^2)*1e-4) +end -function test_twodnavierstokes_deterministicforcingbudgets(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, nμ=0) +function test_twodnavierstokes_deterministicforcing_energybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, nμ=0) n, L = 256, 2π ν, nν = 1e-7, 2 μ, nμ = 1e-1, 0 @@ -120,32 +173,73 @@ function test_twodnavierstokes_deterministicforcingbudgets(dev::Device=CPU(); n= D = Diagnostic(TwoDNavierStokes.energy_dissipation, prob, nsteps=nt) R = Diagnostic(TwoDNavierStokes.energy_drag, prob, nsteps=nt) W = Diagnostic(TwoDNavierStokes.energy_work, prob, nsteps=nt) - Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nt) - DZ = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt) - RZ = Diagnostic(TwoDNavierStokes.enstrophy_drag, prob, nsteps=nt) - WZ = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, nsteps=nt) - diags = [E, D, W, R, Z, DZ, WZ, RZ] + diags = [E, D, W, R] # Step forward stepforward!(prob, diags, nt) TwoDNavierStokes.updatevars!(prob) - E, D, W, R, Z, DZ, WZ, RZ = diags - t = round(μ*cl.t, digits=2) + E, D, W, R = diags i₀ = 1 dEdt = (E[(i₀+1):E.i] - E[i₀:E.i-1])/cl.dt - dZdt = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt ii = (i₀):E.i-1 ii2 = (i₀+1):E.i # dEdt = W - D - R? total = W[ii2] - D[ii] - R[ii] - totalZ = WZ[ii2] - DZ[ii] - RZ[ii] residual = dEdt - total - residualZ = dZdt - totalZ - return isapprox(mean(abs.(residual)), 0, atol=1e-8) && isapprox(mean(abs.(residualZ)), 0, atol=1e-8) + return isapprox(mean(abs.(residual)), 0, atol=1e-8) +end + +function test_twodnavierstokes_deterministicforcing_enstrophybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, nμ=0) + n, L = 256, 2π + ν, nν = 1e-7, 2 + μ, nμ = 1e-1, 0 + dt, tf = 0.005, 0.1/μ + nt = round(Int, tf/dt) + + gr = TwoDGrid(dev, n, L) + x, y = gridpoints(gr) + + # Forcing = 0.01cos(4x)cos(5y)cos(2t) + f = @. 0.01cos(4x) * cos(5y) + fh = rfft(f) + function calcF!(Fh, sol, t, cl, v, p, g::AbstractGrid{T, A}) where {T, A} + Fh = fh * cos(2t) + nothing + end + + prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, nμ=nμ, dt=dt, + stepper="RK4", calcF=calcF!, stochastic=false) + + sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid + + TwoDNavierStokes.set_zeta!(prob, 0*x) + + Z = Diagnostic(TwoDNavierStokes.enstrophy, prob, nsteps=nt) + D = Diagnostic(TwoDNavierStokes.enstrophy_dissipation, prob, nsteps=nt) + R = Diagnostic(TwoDNavierStokes.enstrophy_drag, prob, nsteps=nt) + W = Diagnostic(TwoDNavierStokes.enstrophy_work, prob, nsteps=nt) + diags = [Z, D, W, R] + + # Step forward + stepforward!(prob, diags, nt) + TwoDNavierStokes.updatevars!(prob) + + Z, D, W, R = diags + + i₀ = 1 + dZdt = (Z[(i₀+1):Z.i] - Z[i₀:Z.i-1])/cl.dt + ii = (i₀):Z.i-1 + ii2 = (i₀+1):Z.i + + # dZdt = W - D - R? + total = W[ii2] - D[ii] - R[ii] + residual = dZdt - total + + return isapprox(mean(abs.(residual)), 0, atol=1e-8) end """ From 670c599da4c1827385841f09a51d322a6bd5dcf9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 11 Sep 2020 20:31:54 +1000 Subject: [PATCH 19/19] mentions that enstrophy budgets are computed --- examples/twodnavierstokes_stochasticforcing.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/twodnavierstokes_stochasticforcing.jl b/examples/twodnavierstokes_stochasticforcing.jl index bcb8c455..9468c47b 100644 --- a/examples/twodnavierstokes_stochasticforcing.jl +++ b/examples/twodnavierstokes_stochasticforcing.jl @@ -6,7 +6,8 @@ # A simulation of forced-dissipative two-dimensional turbulence. We solve the # two-dimensional vorticity equation with stochastic excitation and dissipation in # the form of linear drag and hyperviscosity. As a demonstration, we compute how -# each of the forcing and dissipation terms contribute to the energy budget. +# each of the forcing and dissipation terms contribute to the energy and the +# enstrophy budgets. using FourierFlows, Printf, Plots