Functions

Main module for GeophysicalFlows.jl – a collection of solvers for geophysical fluid dynamics problems in periodic domains on CPUs and GPUs. All modules use Fourier-based pseudospectral methods and leverage the functionality of FourierFlows.jl ecosystem.

lambdipole(U, R, grid::TwoDGrid; center=(mean(grid.x), mean(grid.y))

Return the two-dimensional vorticity field of the Lamb dipole with strength U and radius R, centered on center=(xc, yc) and on the grid. Default value for center is the center of the domain.

peakedisotropicspectrum(grid, kpeak, E₀; mask = ones(size(grid.Krsq)), allones = false)

Generate a random two-dimensional relative vorticity field $q(x, y)$ with Fourier spectrum peaked around a central non-dimensional wavenumber kpeak and normalized so that its total kinetic energy is E₀.

Problem(dev::Device = CPU();
                 nx = 256,
                 ny = nx,
                 Lx = 2π,
                 Ly = Lx,
                  ν = 0,
                 nν = 1,
                  μ = 0,
                 nμ = 0,
                 dt = 0.01,
            stepper = "RK4",
              calcF = nothingfunction,
         stochastic = false,
   aliased_fraction = 1/3,
                  T = Float64)

Construct a two-dimensional Navier-Stokes problem on device dev.

Arguments

• dev: (required) CPU() or GPU(); computer architecture used to time-step problem.

Keyword arguments

• nx: Number of grid points in $x$-domain.
• ny: Number of grid points in $y$-domain.
• Lx: Extent of the $x$-domain.
• Ly: Extent of the $y$-domain.
• ν: Small-scale (hyper)-viscosity coefficient.
• nν: (Hyper)-viscosity order, nν$≥ 1$.
• μ: Large-scale (hypo)-viscosity coefficient.
• nμ: (Hypo)-viscosity order, nμ$≤ 0$.
• dt: Time-step.
• stepper: Time-stepping method.
• calcF: Function that calculates the Fourier transform of the forcing, $F̂$.
• stochastic: true or false; boolean denoting whether calcF is temporally stochastic.
• aliased_fraction: the fraction of high-wavenumbers that are zero-ed out by dealias!().
• T: Float32 or Float64; floating point type used for problem data.

energy_dissipation_hyperviscosity(prob)

Return the problem's (prob) domain-averaged energy dissipation rate done by the $ν$ (hyper)-viscosity.

energy_dissipation_hypoviscosity(prob)

Return the problem's (prob) domain-averaged energy dissipation rate done by the $μ$ (hypo)-viscosity.

energy_work(prob)

Return the problem's (prob) domain-averaged rate of work of energy by the forcing $F$,

\[- \int ψ F \frac{𝖽x 𝖽y}{L_x L_y} = - \sum_{𝐤} ψ̂ F̂^* .\]

where $ψ$ is the stream flow.

enstrophy_dissipation_hyperviscosity(prob)

Return the problem's (prob) domain-averaged enstrophy dissipation rate done by the $ν$ (hyper)-viscosity.

enstrophy_dissipation_hypoviscosity(prob)

Return the problem's (prob) domain-averaged enstrophy dissipation rate done by the $μ$ (hypo)-viscosity.

enstrophy_work(prob)

Return the problem's (prob) domain-averaged rate of work of enstrophy by the forcing $F$,

\[\int ζ F \frac{𝖽x 𝖽y}{L_x L_y} = \sum_{𝐤} ζ̂ F̂^* ,\]

where $ζ$ is the relative vorticity.

palinstrophy(prob)

Return the problem's (prob) domain-averaged palinstrophy,

\[\int \frac1{2} |{\bf ∇} ζ|² \frac{𝖽x 𝖽y}{L_x L_y} = \sum_{𝐤} \frac1{2} |𝐤|² |ζ̂|² ,\]

where $ζ$ is the relative vorticity.

calcN_advection!(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term, $- 𝖩(ψ, ζ)$ in conservative form, i.e., $- ∂_x[(∂_y ψ)ζ] - ∂_y[(∂_x ψ)ζ]$ and store it in N:

\[N = - \widehat{𝖩(ψ, ζ)} = - i k_x \widehat{u ζ} - i k_y \widehat{v ζ} .\]

addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term $F̂$ and add it to the nonlinear term $N$.

energy_dissipation(prob, ξ, νξ)

Return the domain-averaged energy dissipation rate done by the viscous term,

\[- ξ (-1)^{n_ξ+1} \int ψ ∇^{2n_ξ} ζ \frac{𝖽x 𝖽y}{L_x L_y} = - ξ \sum_{𝐤} |𝐤|^{2(n_ξ-1)} |ζ̂|² ,\]

where $ξ$ and $nξ$ could be either the (hyper)-viscosity coefficient $ν$ and its order $n_ν$, or the hypo-viscocity coefficient $μ$ and its order $n_μ$.

enstrophy_dissipation(prob, ξ, νξ)

Return the problem's (prob) domain-averaged enstrophy dissipation rate done by the viscous term,

\[ξ (-1)^{n_ξ+1} \int ζ ∇^{2n_ξ} ζ \frac{𝖽x 𝖽y}{L_x L_y} = - ξ \sum_{𝐤} |𝐤|^{2n_ξ} |ζ̂|² ,\]

where $ξ$ and $nξ$ could be either the (hyper)-viscosity coefficient $ν$ and its order $n_ν$, or the hypo-viscocity coefficient $μ$ and its order $n_μ$.

Problem(dev::Device = CPU();
                 nx = 256,
                 ny = nx,
                 Lx = 2π,
                 Ly = Lx,
                  β = 0.0,
 deformation_radius = Inf,
                eta = nothing,
                  ν = 0.0,
                 nν = 1,
                  μ = 0.0,
                 dt = 0.01,
            stepper = "RK4",
              calcF = nothingfunction,
         stochastic = false,
   aliased_fraction = 1/3,
                  T = Float64)

Construct a single-layer quasi-geostrophic problem on device dev.

Arguments

• dev: (required) CPU() or GPU(); computer architecture used to time-step problem.

Keyword arguments

• nx: Number of grid points in $x$-domain.
• ny: Number of grid points in $y$-domain.
• Lx: Extent of the $x$-domain.
• Ly: Extent of the $y$-domain.
• β: Planetary vorticity $y$-gradient.
• deformation_radius: Rossby radius of deformation; set Inf for purely barotropic.
• eta: Topographic potential vorticity.
• ν: Small-scale (hyper)-viscosity coefficient.
• nν: (Hyper)-viscosity order, nν$≥ 1$.
• μ: Linear drag coefficient.
• dt: Time-step.
• stepper: Time-stepping method.
• calcF: Function that calculates the Fourier transform of the forcing, $F̂$.
• stochastic: true or false; boolean denoting whether calcF is temporally stochastic.
• aliased_fraction: the fraction of high-wavenumbers that are zero-ed out by dealias!().
• T: Float32 or Float64; floating point type used for problem data.

streamfunctionfrompv!(ψh, qh, params, grid)

Invert the Fourier transform of PV qh to obtain the Fourier transform of the streamfunction ψh.

energy_dissipation(prob)

Return the problem's (prob) domain-averaged energy dissipation rate.

energy_work(prob)

Return the problem's (prob) domain-averaged rate of work of energy by the forcing F.

energy_drag(prob)

Return the problem's (prob) extraction of domain-averaged energy by drag $μ$.

enstrophy(prob)

Return the problem's (prob) domain-averaged enstrophy

\[\int \frac1{2} (q + η)² \frac{𝖽x 𝖽y}{L_x L_y} = \sum_{𝐤} \frac1{2} |q̂ + η̂|² .\]

enstrophy_dissipation(prob)

Return the problem's (prob) domain-averaged enstrophy dissipation rate.

enstrophy_work(prob)

Return the problem's (prob) domain-averaged rate of work of enstrophy by the forcing $F$.

enstrophy_drag(prob)

Return the problem's (prob) extraction of domain-averaged enstrophy by drag/hypodrag $μ$.

calcN_advection!(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term, $- 𝖩(ψ, q+η)$ in conservative form, i.e., $- ∂_x[(∂_y ψ)(q+η)] - ∂_y[(∂_x ψ)(q+η)]$ and store it in N:

\[N = - \widehat{𝖩(ψ, q + η)} = - i k_x \widehat{u (q + η)} - i k_y \widehat{v (q + η)} .\]

addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term $F̂$ and add it to the nonlinear term $N$.

Problem(nlayers :: Int,
                    dev = CPU();
                     nx = 128,
                     ny = nx,
                     Lx = 2π,
                     Ly = Lx,
                     f₀ = 1.0,
                      β = 0.0,
                      U = zeros(nlayers),
                      H = 1/nlayers * ones(nlayers),
                      b = -(1 .+ 1/nlayers * Array{Float64}(0:nlayers-1)),
                    eta = nothing,
topographic_pv_gradient = (0, 0),
                      μ = 0,
                      ν = 0,
                     nν = 1,
                     dt = 0.01,
                stepper = "RK4",
                 calcFq = nothingfunction,
             stochastic = false,
                 linear = false,
       aliased_fraction = 1/3,
                      T = Float64)

Construct a multi-layer quasi-geostrophic problem with nlayers fluid layers on device dev.

Arguments

• nlayers: (required) Number of fluid layers.
• dev: (required) CPU() (default) or GPU(); computer architecture used to time-step problem.

Keyword arguments

• nx: Number of grid points in $x$-domain.
• ny: Number of grid points in $y$-domain.
• Lx: Extent of the $x$-domain.
• Ly: Extent of the $y$-domain.
• f₀: Constant planetary vorticity.
• β: Planetary vorticity $y$-gradient.
• U: The imposed constant zonal flow $U(y)$ in each fluid layer.
• H: Rest height of each fluid layer.
• b: Boussinesq buoyancy of each fluid layer.
• eta: Periodic component of the topographic potential vorticity.
• topographic_pv_gradient: The $(x, y)$ components of the topographic PV large-scale gradient.
• μ: Linear bottom drag coefficient.
• ν: Small-scale (hyper)-viscosity coefficient.
• nν: (Hyper)-viscosity order, nν$≥ 1$.
• dt: Time-step.
• stepper: Time-stepping method.
• calcF: Function that calculates the Fourier transform of the forcing, $F̂$.
• stochastic: true or false (default); boolean denoting whether calcF is temporally stochastic.
• linear: true or false (default); boolean denoting whether the linearized equations of motions are used.
• aliased_fraction: the fraction of high-wavenumbers that are zero-ed out by dealias!().
• T: Float32 or Float64 (default); floating point type used for problem data.

fwdtransform!(varh, var, params)

Compute the Fourier transform of var and store it in varh.

invtransform!(var, varh, params)

Compute the inverse Fourier transform of varh and store it in var.

streamfunctionfrompv!(ψh, qh, params, grid)

Invert the PV to obtain the Fourier transform of the streamfunction ψh in each layer from qh using ψh = params.S⁻¹ qh.

streamfunctionfrompv!(ψh, qh, params::SingleLayerParams, grid)

Invert the PV to obtain the Fourier transform of the streamfunction ψh for the special case of a single fluid layer configuration. In this case, $ψ̂ = - k⁻² q̂$.

streamfunctionfrompv!(ψh, qh, params::TwoLayerParams, grid)

Invert the PV to obtain the Fourier transform of the streamfunction ψh for the special case of a two fluid layer configuration. In this case we have,

\[ψ̂₁ = - [k² q̂₁ + (f₀² / g′) (q̂₁ / H₂ + q̂₂ / H₁)] / Δ ,\]

\[ψ̂₂ = - [k² q̂₂ + (f₀² / g′) (q̂₁ / H₂ + q̂₂ / H₁)] / Δ ,\]

where $Δ = k² [k² + f₀² (H₁ + H₂) / (g′ H₁ H₂)]$.

(Here, the PV-streamfunction relationship is hard-coded to avoid scalar operations on the GPU.)

pvfromstreamfunction!(qh, ψh, params, grid)

Obtain the Fourier transform of the PV from the streamfunction ψh in each layer using qh = params.S * ψh.

pvfromstreamfunction!(qh, ψh, params::SingleLayerParams, grid)

Obtain the Fourier transform of the PV from the streamfunction ψh for the special case of a single fluid layer configuration. In this case, $q̂ = - k² ψ̂$.

pvfromstreamfunction!(qh, ψh, params::TwoLayerParams, grid)

Obtain the Fourier transform of the PV from the streamfunction ψh for the special case of a two fluid layer configuration. In this case we have,

\[q̂₁ = - k² ψ̂₁ + f₀² / (g′ H₁) (ψ̂₂ - ψ̂₁) ,\]

\[q̂₂ = - k² ψ̂₂ + f₀² / (g′ H₂) (ψ̂₁ - ψ̂₂) .\]

(Here, the PV-streamfunction relationship is hard-coded to avoid scalar operations on the GPU.)

LinearEquation(params, grid)

Return the equation for a multi-layer quasi-geostrophic problem with params and grid. The linear opeartor $L$ includes only (hyper)-viscosity and is computed via hyperviscosity(params, grid).

The nonlinear term is computed via function calcNlinear!.

calcS!(S, Fp, Fm, nlayers, grid)

Construct the array $𝕊$, which consists of nlayer x nlayer static arrays $𝕊_𝐤$ that relate the $q̂_j$'s and $ψ̂_j$'s for every wavenumber: $q̂_𝐤 = 𝕊_𝐤 ψ̂_𝐤$.

calcS⁻¹!(S, Fp, Fm, nlayers, grid)

Construct the array $𝕊⁻¹$, which consists of nlayer x nlayer static arrays $(𝕊_𝐤)⁻¹$ that relate the $q̂_j$'s and $ψ̂_j$'s for every wavenumber: $ψ̂_𝐤 = (𝕊_𝐤)⁻¹ q̂_𝐤$.

calcNlinear!(N, sol, t, clock, vars, params, grid)

Compute the nonlinear term of the linearized equations:

\[N_j = - \widehat{U_j ∂_x Q_j} - \widehat{U_j ∂_x q_j} + \widehat{(∂_y ψ_j)(∂_x Q_j)} - \widehat{(∂_x ψ_j)(∂_y Q_j)} + δ_{j, n} μ |𝐤|^2 ψ̂_n + F̂_j .\]

calcN_advection!(N, sol, vars, params, grid)

Compute the advection term and stores it in N:

\[N_j = - \widehat{𝖩(ψ_j, q_j)} - \widehat{U_j ∂_x Q_j} - \widehat{U_j ∂_x q_j} + \widehat{(∂_y ψ_j)(∂_x Q_j)} - \widehat{(∂_x ψ_j)(∂_y Q_j)} .\]

calcN_linearadvection!(N, sol, vars, params, grid)

Compute the advection term of the linearized equations and stores it in N:

\[N_j = - \widehat{U_j ∂_x Q_j} - \widehat{U_j ∂_x q_j} + \widehat{(∂_y ψ_j)(∂_x Q_j)} - \widehat{(∂_x ψ_j)(∂_y Q_j)} .\]

addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term $F̂$ for each layer and add it to the nonlinear term $N$.

Problem(dev::Device = CPU();
                 nx = 256,
                 Lx = 2π,
                 ny = nx,
                 Ly = Lx,
                  ν = 0,
                 nν = 1,
                 dt = 0.01,
            stepper = "RK4",
              calcF = nothingfunction,
         stochastic = false,
   aliased_fraction = 1/3,
                  T = Float64)

Construct a surface quasi-geostrophic problem on device dev.

Arguments

• dev: (required) CPU() or GPU(); computer architecture used to time-step problem.

Keyword arguments

• nx: Number of grid points in $x$-domain.
• ny: Number of grid points in $y$-domain.
• Lx: Extent of the $x$-domain.
• Ly: Extent of the $y$-domain.
• ν: Small-scale (hyper)-viscosity coefficient.
• nν: (Hyper)-viscosity order, nν$≥ 1$.
• dt: Time-step.
• stepper: Time-stepping method.
• calcF: Function that calculates the Fourier transform of the forcing, $F̂$.
• stochastic: true or false; boolean denoting whether calcF is temporally stochastic.
• aliased_fraction: the fraction of high-wavenumbers that are zero-ed out by dealias!().
• T: Float32 or Float64; floating point type used for problem data.

buoyancy_dissipation(prob)

Return the domain-averaged dissipation rate of surface buoyancy variance due to small scale (hyper)-viscosity,

\[2 ν (-1)^{n_ν} \int b ∇^{2n_ν} b \frac{𝖽x 𝖽y}{L_x L_y} = - 2 ν \sum_{𝐤} |𝐤|^{2n_ν} |b̂|² ,\]

where $ν$ the (hyper)-viscosity coefficient $ν$ and $nν$ the (hyper)-viscosity order. In SQG, this is identical to twice the rate of kinetic energy dissipation

buoyancy_work(prob)
buoyancy_work(sol, vars, grid)

Return the domain-averaged rate of work of buoyancy variance by the forcing,

\[\int 2 b F \frac{𝖽x 𝖽y}{L_x L_y} = \sum_{𝐤} 2 b̂ F̂^* .\]

calcN_advection!(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term, $- 𝖩(ψ, b)$ in conservative form, i.e., $- ∂_x[(∂_y ψ)b] - ∂_y[(∂_x ψ)b]$ and store it in N:

\[N = - \widehat{𝖩(ψ, b)} = - i k_x \widehat{u b} - i k_y \widehat{v b} .\]

addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term $F̂$ and add it to the nonlinear term $N$.

Problem(dev::Device = CPU();
                 nx = 256,
                 ny = nx,
                 Lx = 2π,
                 Ly = Lx,
                  β = 0.0,
                eta = nothing,
                  ν = 0.0,
                 nν = 1,
                  μ = 0.0,
                 dt = 0.01,
            stepper = "RK4",
              calcF = nothingfunction,
         stochastic = false,
   aliased_fraction = 1/3,
                  T = Float64)

Construct a quasi-linear barotropic quasi-geostrophic problem on device dev.

Arguments

• dev: (required) CPU() or GPU(); computer architecture used to time-step problem.

Keyword arguments

• nx: Number of grid points in $x$-domain.
• ny: Number of grid points in $y$-domain.
• Lx: Extent of the $x$-domain.
• Ly: Extent of the $y$-domain.
• β: Planetary vorticity $y$-gradient.
• eta: Topographic potential vorticity.
• ν: Small-scale (hyper)-viscosity coefficient.
• nν: (Hyper)-viscosity order, nν$≥ 1$.
• μ: Linear drag coefficient.
• dt: Time-step.
• stepper: Time-stepping method.
• calcF: Function that calculates the Fourier transform of the forcing, $F̂$.
• stochastic: true or false; boolean denoting whether calcF is temporally stochastic.
• aliased_fraction: the fraction of high-wavenumbers that are zero-ed out by dealias!().
• T: Float32 or Float64; floating point type used for problem data.

dissipation(prob)
dissipation(sol, vars, params, grid)

Return the domain-averaged energy dissipation rate. nν must be >= 1.

work(prob)
work(sol, vars, params, grid)

Return the domain-averaged rate of work of energy by the forcing, params.Fh.

drag(prob)
drag(sol, vars, params, grid)

Return the extraction of domain-averaged energy by drag μ.

calcN_advection!(N, sol, t, clock, vars, params, grid)

Calculate the Fourier transform of the advection term for quasi-linear barotropic QG dynamics.

addforcing!(N, sol, t, clock, vars, params, grid)

When the problem includes forcing, calculate the forcing term $F̂$ and add it to the nonlinear term $N$.