Langmuir turbulence example

This example implements the Langmuir turbulence simulation reported in section 4 of

McWilliams, J. C. et al., "Langmuir Turbulence in the ocean," Journal of Fluid Mechanics (1997).

This example demonstrates:

• how to run large eddy simulations with surface wave effects via the Craik-Leibovich approximation

using Oceananigans

Model set-up

To build the model, we specify the grid, Stokes drift, boundary conditions, and Coriolis parameter.

Domain specification and Grid construction

We create a grid with modest resolution. The grid extent is similar, but not exactly the same as that in McWilliams et al. (1997).

using Oceananigans.Grids

grid = RegularCartesianGrid(size=(32, 32, 48), extent=(128, 128, 96))

RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}
                   domain: x ∈ [0.0, 128.0], y ∈ [0.0, 128.0], z ∈ [-96.0, 0.0]
                 topology: (Periodic, Periodic, Bounded)
  resolution (Nx, Ny, Nz): (32, 32, 48)
   halo size (Hx, Hy, Hz): (1, 1, 1)
grid spacing (Δx, Δy, Δz): (4.0, 4.0, 2.0)

The Stokes Drift profile

The surface wave Stokes drift profile used in McWilliams et al. (1997) corresponds to a 'monochromatic' (that is, single-frequency) wave field with

const wavenumber = 2π / 60 # m⁻¹

and

const amplitude = 0.8 # m

The const declarations ensure that Stokes drift functions compile on the GPU. To run this example on the GPU, write architecture = GPU() in the constructor for IncompressibleModel below.

The Stokes drift at the surface for a monochromatic, deep water wave is

using Oceananigans.Buoyancy: g_Earth

const Uˢ = amplitude^2 * wavenumber * sqrt(g_Earth * wavenumber) # m s⁻¹

0.06791774197745354

The Stokes drift profile is then,

uˢ(z) = Uˢ * exp(2wavenumber * z)

which we need for the initial condition.

Note that Oceananigans.jl implements the Lagrangian-mean form of the Craik-Leibovich equations. This means that our model takes the vertical derivative as an input, rather than the Stokes drift profile itself.

The vertical derivative of the Stokes drift is

∂z_uˢ(z, t) = 2wavenumber * Uˢ * exp(2wavenumber * z)

Finally, we note that the time-derivative of the Stokes drift must be provided if the Stokes drift changes in time. In this example, the Stokes drift is constant and thus the time-derivative of the Stokes drift is 0.

Boundary conditions

At the surface at $z=0$, McWilliams et al. (1997) impose wind stress,

Qᵘ = -3.72e-5 # m² s⁻²

and weak cooling with buoyancy flux

Qᵇ = 2.307e-9 # m³ s⁻²

Oceananigans uses "positive upward" conventions for all fluxes. In consequence, a negative flux at the surface drives positive velocities, and a positive flux of buoyancy drives cooling.

The initial condition and bottom boundary condition for buoyancy impose a linear stratification with buoyancy frequency

N² = 1.936e-5 # s⁻²

To summarize, we impose a surface flux on $u$,

using Oceananigans.BoundaryConditions

u_boundary_conditions = UVelocityBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵘ))

and a surface flux and bottom linear gradient on buoyancy, $b$,

b_boundary_conditions = TracerBoundaryConditions(grid, top = BoundaryCondition(Flux, Qᵇ),
                                                       bottom = BoundaryCondition(Gradient, N²))

Coriolis parameter

McWilliams et al. (1997) use

f = 1e-4 # s⁻¹

which is typical for mid-latitudes on Earth.

Model instantiation

Finally, we are ready to build the model. We use the AnisotropicMinimumDissipation model for large eddy simulation. Because our Stokes drift does not vary in $x, y$, we use UniformStokesDrift, which expects Stokes drift functions of $z, t$ only.

using Oceananigans.Buoyancy: BuoyancyTracer
using Oceananigans.SurfaceWaves: UniformStokesDrift

model = IncompressibleModel(        architecture = CPU(),
                                            grid = grid,
                                         tracers = :b,
                                        buoyancy = BuoyancyTracer(),
                                        coriolis = FPlane(f=f),
                                         closure = AnisotropicMinimumDissipation(),
                                   surface_waves = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                             boundary_conditions = (u=u_boundary_conditions,
                                                    b=b_boundary_conditions),
                            )

IncompressibleModel{CPU, Float64}(time = 0.000 s, iteration = 0) 
├── grid: RegularCartesianGrid{Float64, Periodic, Periodic, Bounded}(Nx=32, Ny=32, Nz=48)
├── tracers: (:b,)
├── closure: VerstappenAnisotropicMinimumDissipation{Float64,NamedTuple{(:b,),Tuple{Float64}},Float64,NamedTuple{(:b,),Tuple{Float64}}}
├── buoyancy: BuoyancyTracer
└── coriolis: FPlane{Float64}

Initial conditions

We make use of random noise concentrated in the upper 4 meters for buoyancy and velocity initial conditions,

Ξ(z) = randn() * exp(z / 4)

Our initial condition for buoyancy consists of a linear stratification, plus noise,

bᵢ(x, y, z) = N² * z + 1e-1 * Ξ(z) * N² * model.grid.Lz

The velocity initial condition is zero Eulerian velocity. This means that we must add the Stokes drift profile to the $u$ velocity field. We also add noise scaled by the friction velocity to $u$ and $w$.

uᵢ(x, y, z) = uˢ(z) + sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

wᵢ(x, y, z) = sqrt(abs(Qᵘ)) * 1e-1 * Ξ(z)

set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

Setting up the simulation

We use the TimeStepWizard for adaptive time-stepping with a Courant-Freidrichs-Lewy (CFL) number of 0.2,

wizard = TimeStepWizard(cfl=0.2, Δt=5.0, max_change=1.1, max_Δt=10.0)

Nice progress messaging

We define a function that prints a helpful message with maximum absolute value of $u, v, w$ and the current wall clock time.

using Oceananigans.Diagnostics, Printf

umax = FieldMaximum(abs, model.velocities.u)
vmax = FieldMaximum(abs, model.velocities.v)
wmax = FieldMaximum(abs, model.velocities.w)

wall_clock = time_ns()

function print_progress(simulation)
    model = simulation.model

    # Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   model.clock.iteration,
                   prettytime(model.clock.time),
                   prettytime(wizard.Δt),
                   umax(), vmax(), wmax(),
                   prettytime(1e-9 * (time_ns() - wall_clock))
                  )

    @info msg

    return nothing
end

print_progress (generic function with 1 method)

Now we create the simulation,

using Oceananigans.Utils: hour # correpsonds to "1 hour", in units of seconds

simulation = Simulation(model, iteration_interval = 100,
                                               Δt = wizard,
                                        stop_time = 4hour,
                                         progress = print_progress)

Simulation{IncompressibleModel{CPU, Float64}}
├── Model clock: time = 0.000 s, iteration = 0 
├── Next time step (TimeStepWizard{Float64}): 5.000 s 
├── Iteration interval: 100
├── Stop criteria: Any[Oceananigans.Simulations.iteration_limit_exceeded, Oceananigans.Simulations.stop_time_exceeded, Oceananigans.Simulations.wall_time_limit_exceeded]
├── Run time: 0.000 s, wall time limit: Inf
├── Stop time: 4.000 hr, stop iteration: Inf
├── Diagnostics: OrderedCollections.OrderedDict with no entries
└── Output writers: OrderedCollections.OrderedDict with no entries

Output

We set up an output writer for the simulation that saves all velocity fields, tracer fields, and the subgrid turbulent diffusivity every 2 minutes.

using Oceananigans.OutputWriters
using Oceananigans.Utils: minute

fields_to_output = merge(model.velocities, model.tracers, (νₑ=model.diffusivities.νₑ,))

simulation.output_writers[:fields] =
    JLD2OutputWriter(model, fields_to_output, time_interval = 2minute,
                     prefix = "langmuir_turbulence", force = true)

Running the simulation

This part is easy,

run!(simulation)

[ Info: i: 0100, t: 8.333 min, Δt: 5.000 s, umax = (8.1e-02, 1.9e-02, 2.6e-02) ms⁻¹, wall time: 17.214 s
[ Info: i: 0200, t: 17.500 min, Δt: 5.500 s, umax = (8.2e-02, 1.9e-02, 2.0e-02) ms⁻¹, wall time: 19.453 s
[ Info: i: 0300, t: 27.583 min, Δt: 6.050 s, umax = (8.5e-02, 2.6e-02, 1.3e-02) ms⁻¹, wall time: 21.687 s
[ Info: i: 0400, t: 38.675 min, Δt: 6.655 s, umax = (8.9e-02, 3.5e-02, 2.0e-02) ms⁻¹, wall time: 23.933 s
[ Info: i: 0500, t: 50.876 min, Δt: 7.321 s, umax = (8.7e-02, 4.8e-02, 2.7e-02) ms⁻¹, wall time: 26.190 s
[ Info: i: 0600, t: 1.072 hr, Δt: 8.053 s, umax = (9.0e-02, 5.2e-02, 2.4e-02) ms⁻¹, wall time: 28.457 s
[ Info: i: 0700, t: 1.318 hr, Δt: 8.858 s, umax = (8.7e-02, 6.9e-02, 2.6e-02) ms⁻¹, wall time: 30.732 s
[ Info: i: 0800, t: 1.573 hr, Δt: 9.207 s, umax = (9.2e-02, 6.3e-02, 3.0e-02) ms⁻¹, wall time: 33.015 s
[ Info: i: 0900, t: 1.815 hr, Δt: 8.692 s, umax = (9.2e-02, 7.1e-02, 2.6e-02) ms⁻¹, wall time: 35.282 s
[ Info: i: 1000, t: 2.055 hr, Δt: 8.658 s, umax = (9.0e-02, 7.1e-02, 3.4e-02) ms⁻¹, wall time: 37.746 s
[ Info: i: 1100, t: 2.301 hr, Δt: 8.860 s, umax = (9.4e-02, 8.0e-02, 3.4e-02) ms⁻¹, wall time: 40.018 s
[ Info: i: 1200, t: 2.538 hr, Δt: 8.507 s, umax = (9.8e-02, 7.3e-02, 2.9e-02) ms⁻¹, wall time: 42.265 s
[ Info: i: 1300, t: 2.764 hr, Δt: 8.133 s, umax = (9.4e-02, 8.2e-02, 2.7e-02) ms⁻¹, wall time: 44.494 s
[ Info: i: 1400, t: 3.000 hr, Δt: 8.504 s, umax = (9.0e-02, 8.3e-02, 2.7e-02) ms⁻¹, wall time: 46.796 s
[ Info: i: 1500, t: 3.246 hr, Δt: 8.866 s, umax = (9.2e-02, 8.9e-02, 2.7e-02) ms⁻¹, wall time: 49.065 s
[ Info: i: 1600, t: 3.489 hr, Δt: 8.737 s, umax = (9.7e-02, 9.2e-02, 3.1e-02) ms⁻¹, wall time: 51.309 s
[ Info: i: 1700, t: 3.719 hr, Δt: 8.277 s, umax = (8.8e-02, 9.6e-02, 3.0e-02) ms⁻¹, wall time: 53.994 s
[ Info: i: 1800, t: 3.951 hr, Δt: 8.372 s, umax = (8.4e-02, 1.0e-01, 3.3e-02) ms⁻¹, wall time: 56.253 s
[ Info: i: 1900, t: 4.167 hr, Δt: 7.755 s, umax = (8.1e-02, 1.0e-01, 3.2e-02) ms⁻¹, wall time: 58.514 s
[ Info: Simulation is stopping. Model time 4.167 hr has hit or exceeded simulation stop time 4.000 hr.

Making a neat movie

We look at the results by plotting vertical slices of $u$ and $w$, and a horizontal slice of $w$ to look for Langmuir cells.

k = searchsortedfirst(grid.zF[:], -8)

Making the coordinate arrays takes a few lines of code,

xw, yw, zw = nodes(model.velocities.w)
xu, yu, zu = nodes(model.velocities.u)

xw, yw, zw = xw[:], yw[:], zw[:]
xu, yu, zu = xu[:], yu[:], zu[:]

Next, we open the JLD2 file, and extract the iterations we ended up saving at,

using JLD2, Plots

file = jldopen(simulation.output_writers[:fields].filepath)

iterations = parse.(Int, keys(file["timeseries/t"]))

WARNING: using Plots.grid in module ex-langmuir_turbulence conflicts with an existing identifier.

This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim)
    levels = range(-clim, stop=clim, length=10)

    cmax = maximum(abs, c)

    if clim < cmax # add levels on either end
        levels = vcat([-cmax], range(-clim, stop=clim, length=10), [cmax])
    end

    return levels
end

Finally, we're ready to animate.

@info "Making an animation from the saved data..."

anim = @animate for (i, iter) in enumerate(iterations)

    @info "Drawing frame $i from iteration $iter \n"

    # Load 3D fields from file
    w = file["timeseries/w/$iter"]
    u = file["timeseries/u/$iter"]

    # Extract slices
    wxy = w[:, :, k]
    wxz = w[:, 1, :]
    uxz = u[:, 1, :]

    wlim = 0.02
    ulim = 0.05
    wlevels = nice_divergent_levels(w, wlim)
    ulevels = nice_divergent_levels(w, ulim)

    wxy_plot = contourf(xw, yw, wxy';
                              color = :balance,
                        aspectratio = :equal,
                              clims = (-wlim, wlim),
                             levels = wlevels,
                              xlims = (0, grid.Lx),
                              ylims = (0, grid.Ly),
                             xlabel = "x (m)",
                             ylabel = "y (m)")

    wxz_plot = contourf(xw, zw, wxz';
                              color = :balance,
                        aspectratio = :equal,
                              clims = (-wlim, wlim),
                             levels = wlevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    uxz_plot = contourf(xu, zu, uxz';
                              color = :balance,
                        aspectratio = :equal,
                              clims = (-ulim, ulim),
                             levels = ulevels,
                              xlims = (0, grid.Lx),
                              ylims = (-grid.Lz, 0),
                             xlabel = "x (m)",
                             ylabel = "z (m)")

    plot(wxy_plot, wxz_plot, uxz_plot, layout=(1, 3), size=(1000, 400),
         title = ["w(x, y, z=-8, t) (m/s)" "w(x, y=0, z, t) (m/s)" "u(x, y = 0, z, t) (m/s)"])

    iter == iterations[end] && close(file)
end

This page was generated using Literate.jl.