From a32f12711ee2b030d75af2c97a0a99c9ae69fd75 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 09:43:50 +1000 Subject: [PATCH 01/43] don't hardcode 3 halo points for pcg --- .../pcg_implicit_free_surface_solver.jl | 35 ++++++++++--------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl b/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl index cfa612d718..9bb5913524 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl @@ -2,7 +2,7 @@ using Oceananigans.Solvers using Oceananigans.Operators using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom using Oceananigans.Architectures -using Oceananigans.Grids: with_halo, isrectilinear +using Oceananigans.Grids: with_halo, isrectilinear, halo_size using Oceananigans.Architectures: device import Oceananigans.Solvers: solve!, precondition! @@ -43,8 +43,8 @@ step `Δt`, gravitational acceleration `g`, and free surface at time-step `n` ` """ function PCGImplicitFreeSurfaceSolver(grid::AbstractGrid, settings, gravitational_acceleration=nothing) # Initialize vertically integrated lateral face areas - ∫ᶻ_Axᶠᶜᶜ = Field((Face, Center, Nothing), with_halo((3, 3, 1), grid)) - ∫ᶻ_Ayᶜᶠᶜ = Field((Center, Face, Nothing), with_halo((3, 3, 1), grid)) + ∫ᶻ_Axᶠᶜᶜ = Field((Face, Center, Nothing), grid) + ∫ᶻ_Ayᶜᶠᶜ = Field((Center, Face, Nothing), grid) vertically_integrated_lateral_areas = (xᶠᶜᶜ = ∫ᶻ_Axᶠᶜᶜ, yᶜᶠᶜ = ∫ᶻ_Ayᶜᶠᶜ) @@ -82,10 +82,10 @@ function solve!(η, implicit_free_surface_solver::PCGImplicitFreeSurfaceSolver, # Take explicit step first? We haven't found improvement from this yet, but perhaps it will # help eventually. #explicit_ab2_step_free_surface!(free_surface, model, Δt, χ) - + ∫ᶻA = implicit_free_surface_solver.vertically_integrated_lateral_areas solver = implicit_free_surface_solver.preconditioned_conjugate_gradient_solver - + # solve!(x, solver, b, args...) solves A*x = b for x. solve!(η, solver, rhs, ∫ᶻA.xᶠᶜᶜ, ∫ᶻA.yᶜᶠᶜ, g, Δt) @@ -118,7 +118,9 @@ end end """ -Returns `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n` + implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) + +Return `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n` and `L` is the linear operator that arises in an implicit time step for the free surface displacement `η`. @@ -142,7 +144,7 @@ end """ _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) -Return the left side of the "implicit η equation" +Return the left side of the "implicit ``η`` equation" ```math (∇ʰ⋅ H ∇ʰ - 1 / (g Δt²)) ηⁿ⁺¹ = 1 / (g Δt) ∇ʰ ⋅ Q★ - 1 / (g Δt²) ηⁿ @@ -157,7 +159,7 @@ and arranged in a symmetric form by multiplying by horizontal areas Az: δⁱÂʷ∂ˣηⁿ⁺¹ + δʲÂˢ∂ʸηⁿ⁺¹ - Az ηⁿ⁺¹ / (g Δt²) = 1 / (g Δt) (δⁱÂʷu̅ˢᵗᵃʳ + δʲÂˢv̅ˢᵗᵃʳ) - Az ηⁿ / (g Δt²) ``` -where ̂ indicates a vertical integral, and +where ̂ indicates a vertical integral, and ̅ indicates a vertical average """ @kernel function _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) @@ -172,20 +174,19 @@ end ##### """ -add to the rhs - H⁻¹ ∇H ⋅ ∇ηⁿ to the rhs... +Add `- H⁻¹ ∇H ⋅ ∇ηⁿ` to the right-hand-side. """ @inline function precondition!(P_r, preconditioner::FFTImplicitFreeSurfaceSolver, r, η, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) poisson_solver = preconditioner.fft_poisson_solver arch = architecture(poisson_solver) grid = preconditioner.three_dimensional_grid Az = grid.Δxᶜᵃᵃ * grid.Δyᵃᶜᵃ # assume horizontal regularity - Lz = grid.Lz + Lz = grid.Lz launch!(arch, grid, :xy, fft_preconditioner_right_hand_side!, poisson_solver.storage, r, η, grid, Az, Lz) - return solve!(P_r, preconditioner, poisson_solver.storage, g, Δt) end @@ -231,20 +232,22 @@ struct DiagonallyDominantInversePreconditioner end """ _diagonally_dominant_precondition!(P_r, grid, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) -Return the diagonally dominant inverse preconditioner applied to the residuals consistently with - `M = D⁻¹(I - (A - D)D⁻¹) ≈ A⁻¹` where `I` is the Identity matrix, -A is the linear operator applied to η, and D is the diagonal of A. +Return the diagonally dominant inverse preconditioner applied to the residuals consistently +with `M = D⁻¹(I - (A - D)D⁻¹) ≈ A⁻¹` where `I` is the identity matrix, `A` is the linear +operator applied to the free surface `η`, and `D` is the diagonal of `A`. ```math P_r = M * r ``` + which expanded in components is + ```math P_rᵢⱼ = rᵢⱼ / Acᵢⱼ - 1 / Acᵢⱼ ( Ax⁻ / Acᵢ₋₁ rᵢ₋₁ⱼ + Ax⁺ / Acᵢ₊₁ rᵢ₊₁ⱼ + Ay⁻ / Acⱼ₋₁ rᵢⱼ₋₁+ Ay⁺ / Acⱼ₊₁ rᵢⱼ₊₁ ) ``` -where `Ac`, `Ax⁻`, `Ax⁺`, `Ay⁻` and `Ay⁺` are the coefficients of -`ηᵢⱼ`, `ηᵢ₋₁ⱼ`, `ηᵢ₊₁ⱼ`, `ηᵢⱼ₋₁` and `ηᵢⱼ₊₁` in `_implicit_free_surface_linear_operation!` +where `Ac`, `Ax⁻`, `Ax⁺`, `Ay⁻` and `Ay⁺` are the coefficients of `ηᵢⱼ`, `ηᵢ₋₁ⱼ`, `ηᵢ₊₁ⱼ`, `ηᵢⱼ₋₁`, +and `ηᵢⱼ₊₁` in `_implicit_free_surface_linear_operation!` """ function diagonally_dominant_precondition!(P_r, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) grid = ∫ᶻ_Axᶠᶜᶜ.grid From e6a12baf5657bce0d722d2af893f972174952f88 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 09:46:22 +1000 Subject: [PATCH 02/43] add barotropic tide example --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index bfa09b4569..b01917d9a0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -31,6 +31,7 @@ examples = [ "One-dimensional diffusion" => "one_dimensional_diffusion", "Two-dimensional turbulence" => "two_dimensional_turbulence", "Internal wave" => "internal_wave", + "Barotropic tide" => "barotropic_tide", "Convecting plankton" => "convecting_plankton", "Ocean wind mixing and convection" => "ocean_wind_mixing_and_convection", "Langmuir turbulence" => "langmuir_turbulence", From d122aaf1d015e4a7da900e13cd20683582c2e866 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 09:46:33 +1000 Subject: [PATCH 03/43] add barotropic tide example --- examples/barotropic_tide.jl | 265 ++++++++++++++++++++++++++++++++++++ 1 file changed, 265 insertions(+) create mode 100644 examples/barotropic_tide.jl diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl new file mode 100644 index 0000000000..6f4e187be0 --- /dev/null +++ b/examples/barotropic_tide.jl @@ -0,0 +1,265 @@ +# # Barotropic tide +# +# In this example, we simulate the evolution of a barotropic tide over a hill. +# +# ## Install dependencies +# +# First let's make sure we have all required packages installed. + +# ```julia +# using Pkg +# pkg"add Oceananigans, CairoMakie" +# ``` + +using Oceananigans +using Oceananigans.Units + +# ## Grid + +# We use an immersed boundary two-dimensional grid (in ``x``-``z``) that is is periodic in +# the ``x``-direction. To construct an immersed boundary grid we first need to create what +# we refer to as "underlying grid", which the grid that encompasses the immersed boundary. + +Nx, Nz = 200, 60 + +H = 2kilometers +Lx = 2000kilometers + +underlying_grid = RectilinearGrid(size = (Nx, Nz), + x = (-Lx/2, Lx/2), + z = (-H, 0), + halo = (4, 4), + topology = (Periodic, Flat, Bounded)) + +# Now we want to create a bathymetry as an immersed boundary grid. + +# ```math +# h(x) = -H + h_0 \\exp(-x^2 / 2σ^2) +# ``` + +h₀ = 50 # m +width = 5kilometers +bump(x, y) = - H + h₀ * exp(-x^2 / 2width^2) + +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bump)) + +# The hill is small; here's how it looks (note that we don't plot all the way to the ocean surface). + +xC = xnodes(grid, Center()) +bottom = grid.immersed_boundary.bottom_height[1:Nx, 1] + +using CairoMakie + +fig = Figure(resolution = (700, 200)) +ax = Axis(fig[1, 1], + xlabel="x [km]", + ylabel="z [m]", + limits=((-Lx/2e3, Lx/2e3), (-H, -4H/5))) + +lines!(ax, xC/1e3, bottom) + +fig + +# Now we want to add a barotropic tide forcing. For example, to add the ``M_2`` tidal forcing +# we need to add forcing in the ``u``-momentum equation of the form: +# ```math +# \\partial_t u = \\dotsb + f_0 \\sin(\\omega_2 t) +# ``` +# where ``\\omega_2 = 2\\pi / T_2``, with ``T_2`` the period of the ``M_2`` tide. + +T₂ = 12.421hours +const ω₂ = 2π / T₂ # radians/sec + +coriolis = FPlane(latitude = -45) + +ε = 0.25 + +U_tidal = ε * ω₂ * width + +const tidal_forcing_amplitude = U_tidal * (coriolis.f^2 - ω₂^2) / ω₂ + +@inline tidal_forcing(x, y, z, t) = tidal_forcing_amplitude * cos(ω₂ * t) + + +# ## Model + +# We built a `HydrostaticFreeSurfaceModel` with an `ImplicitFreeSurface` solver. + +using Oceananigans.Models.HydrostaticFreeSurfaceModels: FFTImplicitFreeSurfaceSolver + +fft_preconditioner = FFTImplicitFreeSurfaceSolver(grid) +free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient, preconditioner=fft_preconditioner) + +model = HydrostaticFreeSurfaceModel(; grid, free_surface, coriolis, + buoyancy = BuoyancyTracer(), + tracers = :b, + momentum_advection = WENO(), + tracer_advection = WENO(), + forcing = (u = tidal_forcing,)) + +# We initialize the model with the tidal flow and a linear stratification. + +uᵢ(x, y, z) = U_tidal + +Nᵢ² = 2e-4 # [s⁻²] initial buoyancy frequency / stratification +bᵢ(x, y, z) = Nᵢ² * z + +set!(model, u=uᵢ, b=bᵢ) + +# Now let's built a `Simulation`. + +Δt = 3minutes +stop_time = 8T₂ + +simulation = Simulation(model, Δt = Δt, stop_time = stop_time) + +# We add a callback to print a message about how the simulation is going, + +using Printf + +wall_clock = Ref(time_ns()) + +function print_progress(sim) + + elapsed = 1e-9 * (time_ns() - wall_clock[]) + + msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹\n", + iteration(sim), prettytime(sim), prettytime(elapsed), + maximum(abs, sim.model.velocities.w)) + + wall_clock[] = time_ns() + + @info msg + + return nothing +end + +simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(200)) + +# ## Diagnostics/Output + +# Add some diagnostics. + +b = model.tracers.b +u, v, w = model.velocities + +U = Field(Average(u)) + +u′ = u - U + +N² = ∂z(b) + +filename = "barotropic_tide" +save_fields_interval = 30minutes + +simulation.output_writers[:fields] = JLD2OutputWriter(model, (; u, u′, w, b, N²); + filename, + schedule = TimeInterval(save_fields_interval), + overwrite_existing = true) + +# We are ready -- let's run! + +@info "Running the simulation..." + +run!(simulation) + +@info "Simulation completed in " * prettytime(simulation.run_wall_time) + +# ## Visualization + +# Now let's visualize our resutls! We use `CairoMakie` in this example. +# On a system with OpenGL `using GLMakie` is more convenient as figures will be +# displayed on the screen. + +using CairoMakie + +# We load the saved buoyancy output on the top, bottom, and east surface as `FieldTimeSeries`es. + +saved_output_filename = filename * ".jld2" + +u′_t = FieldTimeSeries(saved_output_filename, "u′") + w_t = FieldTimeSeries(saved_output_filename, "w") +N²_t = FieldTimeSeries(saved_output_filename, "N²") + +times = u′_t.times +nothing #hide + +# We build the coordinates. We rescale horizontal coordinates so that they correspond to kilometers. + +xu, yu, zu = nodes(u′_t[1]) +xw, yw, zw = nodes(w_t[1]) +xN², yN², zN² = nodes(N²_t[1]) +nothing #hide + +# A utility to mask the region that is within the immersed boundary with `NaN`s. + +using Oceananigans.ImmersedBoundaries: mask_immersed_field! + +function mask_and_get_interior(φ_t, n; value=NaN) + mask_immersed_field!(φ_t[n], value) + return interior(φ_t[n], :, 1, :) +end + +# We use Makie's `Observable` to animate the data. To dive into how `Observable`s work we +# refer to [Makie.jl's Documentation](https://makie.juliaplots.org/stable/documentation/nodes/index.html). + +n = Observable(1) + +title = @lift @sprintf("t = %1.2f days = %1.2f T₂", + round(times[$n] / day, digits=2) , round(times[$n] / T₂, digits=2)) + +u′ₙ = @lift mask_and_get_interior(u′_t, $n) + wₙ = @lift mask_and_get_interior(w_t, $n) +N²ₙ = @lift mask_and_get_interior(N²_t, $n) + +axis_kwargs = (xlabel = "x [km]", + ylabel = "z [m]", + limits = ((-Lx / 2e3, Lx / 2e3), (-H, 0)), + titlesize = 20) + +ulim = 0.8 * maximum(abs, u′_t[end]) +wlim = 0.8 * maximum(abs, w_t[end]) + +fig = Figure(resolution = (700, 900)) + +ax_u = Axis(fig[2, 1]; + title = "u′-velocity", axis_kwargs...) + +ax_w = Axis(fig[3, 1]; + title = "w-velocity", axis_kwargs...) + +ax_N² = Axis(fig[4, 1]; + title = "stratification", axis_kwargs...) + +fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) + +hm_u = heatmap!(ax_u, xu/1e3, zu, u′ₙ; + colorrange = (-ulim, ulim), + colormap = :balance) +Colorbar(fig[2, 2], hm_u) + +hm_w = heatmap!(ax_w, xw/1e3, zw, wₙ; + colorrange = (-wlim, wlim), + colormap = :balance) +Colorbar(fig[3, 2], hm_w) + +hm_N² = heatmap!(ax_N², xN²/1e3, zN², N²ₙ; + colorrange = (0.95Nᵢ², 1.05Nᵢ²), + colormap = :thermal) +Colorbar(fig[4, 2], hm_N²) + +fig + +# Finally, we can record a movie. + +@info "Making an animation from saved data..." + +frames = 1:length(times) + +record(fig, filename * ".mp4", frames, framerate=16) do i + @info string("Plotting frame ", i, " of ", frames[end]) + n[] = i +end +nothing #hide + +# ![](barotropic_tide.mp4) From 292b4e8a3827f938c6d47cb7f1f58c586fe27892 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 09:54:31 +1000 Subject: [PATCH 04/43] run for 4 days --- examples/barotropic_tide.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 6f4e187be0..485e3a8624 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -109,7 +109,7 @@ set!(model, u=uᵢ, b=bᵢ) # Now let's built a `Simulation`. Δt = 3minutes -stop_time = 8T₂ +stop_time = 4days simulation = Simulation(model, Δt = Δt, stop_time = stop_time) @@ -236,17 +236,17 @@ fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) hm_u = heatmap!(ax_u, xu/1e3, zu, u′ₙ; colorrange = (-ulim, ulim), colormap = :balance) -Colorbar(fig[2, 2], hm_u) +Colorbar(fig[2, 2], hm_u, label = "m s⁻¹") hm_w = heatmap!(ax_w, xw/1e3, zw, wₙ; colorrange = (-wlim, wlim), colormap = :balance) -Colorbar(fig[3, 2], hm_w) +Colorbar(fig[3, 2], hm_w, label = "m s⁻¹") hm_N² = heatmap!(ax_N², xN²/1e3, zN², N²ₙ; colorrange = (0.95Nᵢ², 1.05Nᵢ²), colormap = :thermal) -Colorbar(fig[4, 2], hm_N²) +Colorbar(fig[4, 2], hm_N², label = "s⁻²") fig From 72a2399394d60fa22d77dc51822d6175ffcf07e3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 10:00:39 +1000 Subject: [PATCH 05/43] tidal forcing 0 at t=0 --- examples/barotropic_tide.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 485e3a8624..00e24bf4ab 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -78,7 +78,7 @@ U_tidal = ε * ω₂ * width const tidal_forcing_amplitude = U_tidal * (coriolis.f^2 - ω₂^2) / ω₂ -@inline tidal_forcing(x, y, z, t) = tidal_forcing_amplitude * cos(ω₂ * t) +@inline tidal_forcing(x, y, z, t) = tidal_forcing_amplitude * sin(ω₂ * t) # ## Model From 29a44fbc370b8c343a644f07eda78a38cc619c42 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 10:29:51 +1000 Subject: [PATCH 06/43] fix latex rendering --- examples/barotropic_tide.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 00e24bf4ab..91077b474a 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -16,7 +16,7 @@ using Oceananigans.Units # ## Grid -# We use an immersed boundary two-dimensional grid (in ``x``-``z``) that is is periodic in +# We use an immersed boundary two-dimensional grid (in ``x``--``z``) that is is periodic in # the ``x``-direction. To construct an immersed boundary grid we first need to create what # we refer to as "underlying grid", which the grid that encompasses the immersed boundary. @@ -34,7 +34,7 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), # Now we want to create a bathymetry as an immersed boundary grid. # ```math -# h(x) = -H + h_0 \\exp(-x^2 / 2σ^2) +# h(x) = -H + h_0 \exp(-x^2 / 2σ^2) # ``` h₀ = 50 # m @@ -63,9 +63,9 @@ fig # Now we want to add a barotropic tide forcing. For example, to add the ``M_2`` tidal forcing # we need to add forcing in the ``u``-momentum equation of the form: # ```math -# \\partial_t u = \\dotsb + f_0 \\sin(\\omega_2 t) +# \partial_t u = \dotsb + f_0 \sin(\omega_2 t) # ``` -# where ``\\omega_2 = 2\\pi / T_2``, with ``T_2`` the period of the ``M_2`` tide. +# where ``\omega_2 = 2π / T_2``, with ``T_2`` the period of the ``M_2`` tide. T₂ = 12.421hours const ω₂ = 2π / T₂ # radians/sec From d2258d8770d966fcf9a6c94e1056073bd5d6b757 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 10:34:32 +1000 Subject: [PATCH 07/43] adds remarks about saved output --- examples/barotropic_tide.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 91077b474a..56f8b09b97 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -138,7 +138,9 @@ simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterv # ## Diagnostics/Output -# Add some diagnostics. +# Add some diagnostics. Instead of ``u`` we save the deviation of ``u`` from its instantaneous +# domain average, ``u' = u - (L_x H)^{-1} \int u \, \mathrm{d}x \mathrm{d}z``. We also save +# the stratification ``N^2 = \partial_z b``. b = model.tracers.b u, v, w = model.velocities @@ -173,7 +175,7 @@ run!(simulation) using CairoMakie -# We load the saved buoyancy output on the top, bottom, and east surface as `FieldTimeSeries`es. +# We load the saved valocity and stratification output as `FieldTimeSeries`es. saved_output_filename = filename * ".jld2" From bbd8ee43b96506ac437647e66974f12c4bcd6597 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:00:17 +1000 Subject: [PATCH 08/43] add more explanation regarding the tidal forcing --- examples/barotropic_tide.jl | 36 +++++++++++++++++++++++++++++------- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 56f8b09b97..71f5f269c8 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -63,24 +63,46 @@ fig # Now we want to add a barotropic tide forcing. For example, to add the ``M_2`` tidal forcing # we need to add forcing in the ``u``-momentum equation of the form: # ```math -# \partial_t u = \dotsb + f_0 \sin(\omega_2 t) +# \partial_t u = \dotsb + F_0 \sin(\omega_2 t) # ``` -# where ``\omega_2 = 2π / T_2``, with ``T_2`` the period of the ``M_2`` tide. +# where ``\omega_2 = 2π / T_2``, with ``T_2 = 12.421 \,\mathrm{hours}`` the period of the ``M_2`` tide. -T₂ = 12.421hours -const ω₂ = 2π / T₂ # radians/sec +# The excursion parameter is a nondimensional number that expresses the ratio of the flow movement +# due to the tide compared to the size of the width of the hill. +# +# ```math +# \epsilon = \frac{U_{\mathrm{tidal}} / \omega_2}{\sigma} +# ``` +# +# We prescribe the excursion parameter which, in turn, implies a tidal velocity ``U_{\mathrm{tidal}}`` +# which then allows us to determing the tidal forcing amplitude ``F_0``. For the last step, we +# use Fourier decomposition on the inviscid, linearized momentum equations to determine the +# flow response for a given tidal forcing. Doing so we get that for the sinusoidal forcing above, +# the tidal velocity and tidal forcing amplitudes are related via: +# +# ```math +# U_{\mathrm{tidal}} = \frac{\omega F_0}{\omega^2 - f^2} +# ``` +# +# The Coriolis frequency is needed, so we start by constructing a Coriolis on an ``f``-plane at the +# mid-latitudes. coriolis = FPlane(latitude = -45) -ε = 0.25 +# Now we have everything we require to construct the tidal forcing given a value of the +# excursion parameter. +# +T₂ = 12.421hours +const ω₂ = 2π / T₂ # radians/sec + +ε = 0.25 # the excursion parameter U_tidal = ε * ω₂ * width -const tidal_forcing_amplitude = U_tidal * (coriolis.f^2 - ω₂^2) / ω₂ +const tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂ @inline tidal_forcing(x, y, z, t) = tidal_forcing_amplitude * sin(ω₂ * t) - # ## Model # We built a `HydrostaticFreeSurfaceModel` with an `ImplicitFreeSurface` solver. From 900db08dd80201988949ea7aa1bf93d63bbc74b3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:09:12 +1000 Subject: [PATCH 09/43] add remarks about IBG --- examples/barotropic_tide.jl | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 71f5f269c8..10ec03dc58 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -31,7 +31,11 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4), topology = (Periodic, Flat, Bounded)) -# Now we want to create a bathymetry as an immersed boundary grid. +# Now we can create the non-trivial bathymetry as an immersed boundary. We use `GridFittedBottom` that +# gets as input a two-dimensional function whose arguments are the grid's native horizontal coordinates, +# that returns the depth of the bottom. +# +# In this example we'd like to have a small Gaussian hill at the center of the domain. # ```math # h(x) = -H + h_0 \exp(-x^2 / 2σ^2) @@ -39,9 +43,10 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), h₀ = 50 # m width = 5kilometers -bump(x, y) = - H + h₀ * exp(-x^2 / 2width^2) +hill(x) = h₀ * exp(-x^2 / 2width^2) +bottom(x, y) = - H + hill(x) -grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bump)) +grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom)) # The hill is small; here's how it looks (note that we don't plot all the way to the ocean surface). From e067cea4623847bbc18b9386d9475fb4bfc98a44 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:11:39 +1000 Subject: [PATCH 10/43] more remarks about IBG --- examples/barotropic_tide.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 10ec03dc58..6973aaa539 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -32,8 +32,9 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), topology = (Periodic, Flat, Bounded)) # Now we can create the non-trivial bathymetry as an immersed boundary. We use `GridFittedBottom` that -# gets as input a two-dimensional function whose arguments are the grid's native horizontal coordinates, -# that returns the depth of the bottom. +# gets as input either *(i)* a two-dimensional function whose arguments are the grid's native horizontal +# coordinates and it returns the ``z`` of the bottom, or *(i)* a two-dimensional array with the values +# of ``z`` at the bottom at the cell centers. # # In this example we'd like to have a small Gaussian hill at the center of the domain. From 8cbe58877ff70a6dfbfa44af23acf6864afed776 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:12:14 +1000 Subject: [PATCH 11/43] (i) -> (ii) --- examples/barotropic_tide.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 6973aaa539..82d211c8b0 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -33,7 +33,7 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), # Now we can create the non-trivial bathymetry as an immersed boundary. We use `GridFittedBottom` that # gets as input either *(i)* a two-dimensional function whose arguments are the grid's native horizontal -# coordinates and it returns the ``z`` of the bottom, or *(i)* a two-dimensional array with the values +# coordinates and it returns the ``z`` of the bottom, or *(ii)* a two-dimensional array with the values # of ``z`` at the bottom at the cell centers. # # In this example we'd like to have a small Gaussian hill at the center of the domain. From ba801a2b13e1b120149128f210167830c7177db8 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:15:46 +1000 Subject: [PATCH 12/43] first load, then plot --- examples/barotropic_tide.jl | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 82d211c8b0..85922a978c 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -195,15 +195,9 @@ run!(simulation) @info "Simulation completed in " * prettytime(simulation.run_wall_time) -# ## Visualization +# ## Load output -# Now let's visualize our resutls! We use `CairoMakie` in this example. -# On a system with OpenGL `using GLMakie` is more convenient as figures will be -# displayed on the screen. - -using CairoMakie - -# We load the saved valocity and stratification output as `FieldTimeSeries`es. +# First, we load the saved valocity and stratification output as `FieldTimeSeries`es. saved_output_filename = filename * ".jld2" @@ -214,14 +208,22 @@ N²_t = FieldTimeSeries(saved_output_filename, "N²") times = u′_t.times nothing #hide -# We build the coordinates. We rescale horizontal coordinates so that they correspond to kilometers. +# We retrieve each field's coordinates. xu, yu, zu = nodes(u′_t[1]) xw, yw, zw = nodes(w_t[1]) xN², yN², zN² = nodes(N²_t[1]) nothing #hide -# A utility to mask the region that is within the immersed boundary with `NaN`s. +# ## Visualize + +# Now we can visualize our resutls! We use `CairoMakie` here. On a system with OpenGL +# `using GLMakie` is more convenient as figures will be displayed on the screen. + +using CairoMakie + +# First, a utility to mask the region that is within the immersed boundary with `NaN`s and +# also select the ``x``-``z`` interior slices of our fields. using Oceananigans.ImmersedBoundaries: mask_immersed_field! From b3f39a47eaf37b23e5901f1159a28788d9f85d0c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:20:31 +1000 Subject: [PATCH 13/43] missing _2 --- examples/barotropic_tide.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 85922a978c..026178d89c 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -87,7 +87,7 @@ fig # the tidal velocity and tidal forcing amplitudes are related via: # # ```math -# U_{\mathrm{tidal}} = \frac{\omega F_0}{\omega^2 - f^2} +# U_{\mathrm{tidal}} = \frac{\omega_2 F_0}{\omega_2^2 - f^2} # ``` # # The Coriolis frequency is needed, so we start by constructing a Coriolis on an ``f``-plane at the From 803d21bb9efb8f865b10212fc07a1b45d7826dc1 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:22:17 +1000 Subject: [PATCH 14/43] some more remarks --- examples/barotropic_tide.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 026178d89c..50d351ba42 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -87,10 +87,11 @@ fig # the tidal velocity and tidal forcing amplitudes are related via: # # ```math -# U_{\mathrm{tidal}} = \frac{\omega_2 F_0}{\omega_2^2 - f^2} +# U_{\mathrm{tidal}} = \frac{\omega_2}{\omega_2^2 - f^2}F_0 # ``` # -# The Coriolis frequency is needed, so we start by constructing a Coriolis on an ``f``-plane at the +# Now we have the way to find the value of the tidal forcing amplitude that would correspond to a +# given excursion parameter. The Coriolis frequency is needed, so we start by constructing a Coriolis on an ``f``-plane at the # mid-latitudes. coriolis = FPlane(latitude = -45) From dc0ea6336b181c97ed468fe1a314023f775b023c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:28:37 +1000 Subject: [PATCH 15/43] fix bottom redefinition + beautifications --- examples/barotropic_tide.jl | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 50d351ba42..e50f0f3282 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -20,7 +20,7 @@ using Oceananigans.Units # the ``x``-direction. To construct an immersed boundary grid we first need to create what # we refer to as "underlying grid", which the grid that encompasses the immersed boundary. -Nx, Nz = 200, 60 +Nx, Nz = 200, 80 H = 2kilometers Lx = 2000kilometers @@ -52,7 +52,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom)) # The hill is small; here's how it looks (note that we don't plot all the way to the ocean surface). xC = xnodes(grid, Center()) -bottom = grid.immersed_boundary.bottom_height[1:Nx, 1] +bottom_boundary = grid.immersed_boundary.bottom_height[1:Nx, 1] using CairoMakie @@ -62,7 +62,7 @@ ax = Axis(fig[1, 1], ylabel="z [m]", limits=((-Lx/2e3, Lx/2e3), (-H, -4H/5))) -lines!(ax, xC/1e3, bottom) +lines!(ax, xC/1e3, bottom_boundary) fig @@ -109,6 +109,7 @@ U_tidal = ε * ω₂ * width const tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂ @inline tidal_forcing(x, y, z, t) = tidal_forcing_amplitude * sin(ω₂ * t) +nothing #hide # ## Model @@ -164,6 +165,7 @@ function print_progress(sim) end simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(200)) +nothing #hide # ## Diagnostics/Output @@ -190,12 +192,8 @@ simulation.output_writers[:fields] = JLD2OutputWriter(model, (; u, u′, w, b, N # We are ready -- let's run! -@info "Running the simulation..." - run!(simulation) -@info "Simulation completed in " * prettytime(simulation.run_wall_time) - # ## Load output # First, we load the saved valocity and stratification output as `FieldTimeSeries`es. @@ -232,6 +230,7 @@ function mask_and_get_interior(φ_t, n; value=NaN) mask_immersed_field!(φ_t[n], value) return interior(φ_t[n], :, 1, :) end +nothing #hide # We use Makie's `Observable` to animate the data. To dive into how `Observable`s work we # refer to [Makie.jl's Documentation](https://makie.juliaplots.org/stable/documentation/nodes/index.html). @@ -256,13 +255,13 @@ wlim = 0.8 * maximum(abs, w_t[end]) fig = Figure(resolution = (700, 900)) ax_u = Axis(fig[2, 1]; - title = "u′-velocity", axis_kwargs...) + title = L"$u'$-velocity", axis_kwargs...) ax_w = Axis(fig[3, 1]; - title = "w-velocity", axis_kwargs...) + title = L"$w$-velocity", axis_kwargs...) ax_N² = Axis(fig[4, 1]; - title = "stratification", axis_kwargs...) + title = L"stratification $N^2$", axis_kwargs...) fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) From fbde03b2d8bc29d9d56e052876053507285c9cff Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:33:44 +1000 Subject: [PATCH 16/43] simpler panel titles --- examples/barotropic_tide.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index e50f0f3282..fe59ce98c7 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -255,13 +255,13 @@ wlim = 0.8 * maximum(abs, w_t[end]) fig = Figure(resolution = (700, 900)) ax_u = Axis(fig[2, 1]; - title = L"$u'$-velocity", axis_kwargs...) + title = "u'-velocity", axis_kwargs...) ax_w = Axis(fig[3, 1]; - title = L"$w$-velocity", axis_kwargs...) + title = "w-velocity", axis_kwargs...) ax_N² = Axis(fig[4, 1]; - title = L"stratification $N^2$", axis_kwargs...) + title = "stratification N²", axis_kwargs...) fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) From 373e5f518f4c79cf82f77d06556075bd9b788fc6 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:38:26 +1000 Subject: [PATCH 17/43] simpler plotting --- examples/barotropic_tide.jl | 34 +++++++++++----------------------- 1 file changed, 11 insertions(+), 23 deletions(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index fe59ce98c7..806c03cf99 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -207,11 +207,11 @@ N²_t = FieldTimeSeries(saved_output_filename, "N²") times = u′_t.times nothing #hide -# We retrieve each field's coordinates. +# We retrieve each field's coordinates and convert from meters to kilometers. -xu, yu, zu = nodes(u′_t[1]) -xw, yw, zw = nodes(w_t[1]) -xN², yN², zN² = nodes(N²_t[1]) +xu, yu, zu = nodes(u′_t[1]) ./ 1e3 +xw, yw, zw = nodes(w_t[1]) ./ 1e3 +xN², yN², zN² = nodes(N²_t[1]) ./ 1e3 nothing #hide # ## Visualize @@ -246,7 +246,7 @@ N²ₙ = @lift mask_and_get_interior(N²_t, $n) axis_kwargs = (xlabel = "x [km]", ylabel = "z [m]", - limits = ((-Lx / 2e3, Lx / 2e3), (-H, 0)), + limits = ((-Lx / 2e3, Lx / 2e3), (-H / 1e3, 0)), # note conversion to kilometers titlesize = 20) ulim = 0.8 * maximum(abs, u′_t[end]) @@ -254,30 +254,18 @@ wlim = 0.8 * maximum(abs, w_t[end]) fig = Figure(resolution = (700, 900)) -ax_u = Axis(fig[2, 1]; - title = "u'-velocity", axis_kwargs...) - -ax_w = Axis(fig[3, 1]; - title = "w-velocity", axis_kwargs...) - -ax_N² = Axis(fig[4, 1]; - title = "stratification N²", axis_kwargs...) - fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) -hm_u = heatmap!(ax_u, xu/1e3, zu, u′ₙ; - colorrange = (-ulim, ulim), - colormap = :balance) +ax_u = Axis(fig[2, 1]; title = "u'-velocity", axis_kwargs...) +hm_u = heatmap!(ax_u, xu, zu, u′ₙ; colorrange = (-ulim, ulim), colormap = :balance) Colorbar(fig[2, 2], hm_u, label = "m s⁻¹") -hm_w = heatmap!(ax_w, xw/1e3, zw, wₙ; - colorrange = (-wlim, wlim), - colormap = :balance) +ax_w = Axis(fig[3, 1]; title = "w-velocity", axis_kwargs...) +hm_w = heatmap!(ax_w, xw, zw, wₙ; colorrange = (-wlim, wlim), colormap = :balance) Colorbar(fig[3, 2], hm_w, label = "m s⁻¹") -hm_N² = heatmap!(ax_N², xN²/1e3, zN², N²ₙ; - colorrange = (0.95Nᵢ², 1.05Nᵢ²), - colormap = :thermal) +ax_N² = Axis(fig[4, 1]; title = "stratification N²", axis_kwargs...) +hm_N² = heatmap!(ax_N², xN², zN², N²ₙ; colorrange = (0.95Nᵢ², 1.05Nᵢ²), colormap = :thermal) Colorbar(fig[4, 2], hm_N², label = "s⁻²") fig From 9f6590dcc861bf90b24b3779bc322d598d083d70 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:39:22 +1000 Subject: [PATCH 18/43] z [m] -> km --- examples/barotropic_tide.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index 806c03cf99..c6363652c4 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -245,7 +245,7 @@ u′ₙ = @lift mask_and_get_interior(u′_t, $n) N²ₙ = @lift mask_and_get_interior(N²_t, $n) axis_kwargs = (xlabel = "x [km]", - ylabel = "z [m]", + ylabel = "z [km]", limits = ((-Lx / 2e3, Lx / 2e3), (-H / 1e3, 0)), # note conversion to kilometers titlesize = 20) From b8c0efb54d699994bc4286b09c48acdd6afe85bb Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 4 Jun 2023 11:44:37 +1000 Subject: [PATCH 19/43] simplify --- examples/barotropic_tide.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/barotropic_tide.jl b/examples/barotropic_tide.jl index c6363652c4..c6095721e1 100644 --- a/examples/barotropic_tide.jl +++ b/examples/barotropic_tide.jl @@ -141,7 +141,7 @@ set!(model, u=uᵢ, b=bᵢ) Δt = 3minutes stop_time = 4days -simulation = Simulation(model, Δt = Δt, stop_time = stop_time) +simulation = Simulation(model; Δt, stop_time) # We add a callback to print a message about how the simulation is going, From 86d622f1ad9d5ff861cc5c6792f6937c1a831877 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 5 Jun 2023 13:53:35 +1000 Subject: [PATCH 20/43] rename to internal tide + add some remarks from review --- docs/make.jl | 30 ++++++++-------- .../{barotropic_tide.jl => internal_tide.jl} | 35 +++++++++---------- 2 files changed, 32 insertions(+), 33 deletions(-) rename examples/{barotropic_tide.jl => internal_tide.jl} (86%) diff --git a/docs/make.jl b/docs/make.jl index b01917d9a0..3b2d62160c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,18 +28,18 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/generated") examples = [ - "One-dimensional diffusion" => "one_dimensional_diffusion", - "Two-dimensional turbulence" => "two_dimensional_turbulence", - "Internal wave" => "internal_wave", - "Barotropic tide" => "barotropic_tide", - "Convecting plankton" => "convecting_plankton", - "Ocean wind mixing and convection" => "ocean_wind_mixing_and_convection", - "Langmuir turbulence" => "langmuir_turbulence", - "Baroclinic adjustment" => "baroclinic_adjustment", - "Kelvin-Helmholtz instability" => "kelvin_helmholtz_instability", - "Shallow water Bickley jet" => "shallow_water_Bickley_jet", - "Horizontal convection" => "horizontal_convection", - "Tilted bottom boundary layer" => "tilted_bottom_boundary_layer" + # "One-dimensional diffusion" => "one_dimensional_diffusion", + # "Two-dimensional turbulence" => "two_dimensional_turbulence", + # "Internal wave" => "internal_wave", + "Internal tide by a seamount" => "internal_tide", + # "Convecting plankton" => "convecting_plankton", + # "Ocean wind mixing and convection" => "ocean_wind_mixing_and_convection", + # "Langmuir turbulence" => "langmuir_turbulence", + # "Baroclinic adjustment" => "baroclinic_adjustment", + # "Kelvin-Helmholtz instability" => "kelvin_helmholtz_instability", + # "Shallow water Bickley jet" => "shallow_water_Bickley_jet", + # "Horizontal convection" => "horizontal_convection", + # "Tilted bottom boundary layer" => "tilted_bottom_boundary_layer" ] example_scripts = [ filename * ".jl" for (title, filename) in examples ] @@ -144,10 +144,10 @@ makedocs(bib, format = format, pages = pages, modules = [Oceananigans], - doctest = true, - strict = true, + doctest = false, + strict = false, clean = true, - checkdocs = :exports + checkdocs = :none ) @info "Clean up temporary .jld2 and .nc output created by doctests or literated examples..." diff --git a/examples/barotropic_tide.jl b/examples/internal_tide.jl similarity index 86% rename from examples/barotropic_tide.jl rename to examples/internal_tide.jl index c6095721e1..5a485c6b32 100644 --- a/examples/barotropic_tide.jl +++ b/examples/internal_tide.jl @@ -16,9 +16,8 @@ using Oceananigans.Units # ## Grid -# We use an immersed boundary two-dimensional grid (in ``x``--``z``) that is is periodic in -# the ``x``-direction. To construct an immersed boundary grid we first need to create what -# we refer to as "underlying grid", which the grid that encompasses the immersed boundary. +# We create an `ImmersedBoundaryGrid` wrapped around an underlying two-dimensional `RectilinearGrid` +# that is periodic in ``x`` and bounded in ``z``. Nx, Nz = 200, 80 @@ -31,13 +30,13 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), halo = (4, 4), topology = (Periodic, Flat, Bounded)) -# Now we can create the non-trivial bathymetry as an immersed boundary. We use `GridFittedBottom` that -# gets as input either *(i)* a two-dimensional function whose arguments are the grid's native horizontal -# coordinates and it returns the ``z`` of the bottom, or *(ii)* a two-dimensional array with the values -# of ``z`` at the bottom at the cell centers. +# Now we can create the non-trivial bathymetry. We use `GridFittedBottom` that gets as input either +# *(i)* a two-dimensional function whose arguments are the grid's native horizontal coordinates and +# it returns the ``z`` of the bottom, or *(ii)* a two-dimensional array with the values of ``z`` at +# the bottom cell centers. # # In this example we'd like to have a small Gaussian hill at the center of the domain. - +# # ```math # h(x) = -H + h_0 \exp(-x^2 / 2σ^2) # ``` @@ -113,14 +112,14 @@ nothing #hide # ## Model -# We built a `HydrostaticFreeSurfaceModel` with an `ImplicitFreeSurface` solver. - -using Oceananigans.Models.HydrostaticFreeSurfaceModels: FFTImplicitFreeSurfaceSolver +# We built a `HydrostaticFreeSurfaceModel` with a `SplitExplicitFreeSurface` solver. +# We limit our maximum timestep to 3 minutes. -fft_preconditioner = FFTImplicitFreeSurfaceSolver(grid) -free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient, preconditioner=fft_preconditioner) +max_Δt = 3minutes +free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7, max_Δt) -model = HydrostaticFreeSurfaceModel(; grid, free_surface, coriolis, +model = HydrostaticFreeSurfaceModel(; grid, coriolis, + free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7, max_Δt), buoyancy = BuoyancyTracer(), tracers = :b, momentum_advection = WENO(), @@ -138,7 +137,7 @@ set!(model, u=uᵢ, b=bᵢ) # Now let's built a `Simulation`. -Δt = 3minutes +Δt = max_Δt stop_time = 4days simulation = Simulation(model; Δt, stop_time) @@ -182,7 +181,7 @@ u′ = u - U N² = ∂z(b) -filename = "barotropic_tide" +filename = "internal_tide" save_fields_interval = 30minutes simulation.output_writers[:fields] = JLD2OutputWriter(model, (; u, u′, w, b, N²); @@ -196,7 +195,7 @@ run!(simulation) # ## Load output -# First, we load the saved valocity and stratification output as `FieldTimeSeries`es. +# First, we load the saved velocities and stratification output as `FieldTimeSeries`es. saved_output_filename = filename * ".jld2" @@ -282,4 +281,4 @@ record(fig, filename * ".mp4", frames, framerate=16) do i end nothing #hide -# ![](barotropic_tide.mp4) +# ![](internal_tide.mp4) From 2f8c059dafc30702c599d0e046caddd75badeca4 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 5 Jun 2023 22:01:35 +1000 Subject: [PATCH 21/43] apply review remarks --- examples/internal_tide.jl | 71 ++++++++++++++++----------------------- 1 file changed, 29 insertions(+), 42 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 5a485c6b32..6b6b71a752 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -1,6 +1,7 @@ -# # Barotropic tide +# # Internal tide by a seamount # -# In this example, we simulate the evolution of a barotropic tide over a hill. +# In this example, we show how internal tide is generated from a barotropic tidal flow +# sloshing back and forth over a sea mount. # # ## Install dependencies # @@ -19,14 +20,11 @@ using Oceananigans.Units # We create an `ImmersedBoundaryGrid` wrapped around an underlying two-dimensional `RectilinearGrid` # that is periodic in ``x`` and bounded in ``z``. -Nx, Nz = 200, 80 - -H = 2kilometers -Lx = 2000kilometers +Nx, Nz = 250, 125 underlying_grid = RectilinearGrid(size = (Nx, Nz), - x = (-Lx/2, Lx/2), - z = (-H, 0), + x = (-1000kilometers, 1000kilometers), + z = (-2kilometers, 0), halo = (4, 4), topology = (Periodic, Flat, Bounded)) @@ -41,8 +39,8 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), # h(x) = -H + h_0 \exp(-x^2 / 2σ^2) # ``` -h₀ = 50 # m -width = 5kilometers +h₀ = 250 # m +width = 20kilometers hill(x) = h₀ * exp(-x^2 / 2width^2) bottom(x, y) = - H + hill(x) @@ -59,7 +57,7 @@ fig = Figure(resolution = (700, 200)) ax = Axis(fig[1, 1], xlabel="x [km]", ylabel="z [m]", - limits=((-Lx/2e3, Lx/2e3), (-H, -4H/5))) + limits=((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz, -4grid.Lz/5))) lines!(ax, xC/1e3, bottom_boundary) @@ -99,38 +97,38 @@ coriolis = FPlane(latitude = -45) # excursion parameter. # T₂ = 12.421hours -const ω₂ = 2π / T₂ # radians/sec +ω₂ = 2π / T₂ # radians/sec -ε = 0.25 # the excursion parameter +ε = 0.1 # excursion parameter U_tidal = ε * ω₂ * width -const tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂ +tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂ -@inline tidal_forcing(x, y, z, t) = tidal_forcing_amplitude * sin(ω₂ * t) -nothing #hide +@inline tidal_forcing(x, y, z, t, p) = p.tidal_forcing_amplitude * sin(p.ω₂ * t) + +u_forcing = Forcing(tidal_forcing, parameters=(; tidal_forcing_amplitude, ω₂)) # ## Model # We built a `HydrostaticFreeSurfaceModel` with a `SplitExplicitFreeSurface` solver. # We limit our maximum timestep to 3 minutes. -max_Δt = 3minutes +max_Δt = 5minutes free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7, max_Δt) -model = HydrostaticFreeSurfaceModel(; grid, coriolis, - free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7, max_Δt), +model = HydrostaticFreeSurfaceModel(; grid, coriolis, free_surface, buoyancy = BuoyancyTracer(), tracers = :b, momentum_advection = WENO(), tracer_advection = WENO(), - forcing = (u = tidal_forcing,)) + forcing = (; u = u_forcing)) # We initialize the model with the tidal flow and a linear stratification. uᵢ(x, y, z) = U_tidal -Nᵢ² = 2e-4 # [s⁻²] initial buoyancy frequency / stratification +Nᵢ² = 1e-4 # [s⁻²] initial buoyancy frequency / stratification bᵢ(x, y, z) = Nᵢ² * z set!(model, u=uᵢ, b=bᵢ) @@ -217,39 +215,28 @@ nothing #hide # Now we can visualize our resutls! We use `CairoMakie` here. On a system with OpenGL # `using GLMakie` is more convenient as figures will be displayed on the screen. - -using CairoMakie - -# First, a utility to mask the region that is within the immersed boundary with `NaN`s and -# also select the ``x``-``z`` interior slices of our fields. - -using Oceananigans.ImmersedBoundaries: mask_immersed_field! - -function mask_and_get_interior(φ_t, n; value=NaN) - mask_immersed_field!(φ_t[n], value) - return interior(φ_t[n], :, 1, :) -end -nothing #hide - +# # We use Makie's `Observable` to animate the data. To dive into how `Observable`s work we # refer to [Makie.jl's Documentation](https://makie.juliaplots.org/stable/documentation/nodes/index.html). +using CairoMakie + n = Observable(1) title = @lift @sprintf("t = %1.2f days = %1.2f T₂", round(times[$n] / day, digits=2) , round(times[$n] / T₂, digits=2)) -u′ₙ = @lift mask_and_get_interior(u′_t, $n) - wₙ = @lift mask_and_get_interior(w_t, $n) -N²ₙ = @lift mask_and_get_interior(N²_t, $n) +u′ₙ = @lift interior(u′_t[$n], :, 1, :) + wₙ = @lift interior( w_t[$n], :, 1, :) +N²ₙ = @lift interior(N²_t[$n], :, 1, :) axis_kwargs = (xlabel = "x [km]", ylabel = "z [km]", - limits = ((-Lx / 2e3, Lx / 2e3), (-H / 1e3, 0)), # note conversion to kilometers + limits = ((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz/1e3, 0)), # note conversion to kilometers titlesize = 20) -ulim = 0.8 * maximum(abs, u′_t[end]) -wlim = 0.8 * maximum(abs, w_t[end]) +ulim = maximum(abs, u′_t[end]) +wlim = maximum(abs, w_t[end]) fig = Figure(resolution = (700, 900)) @@ -264,7 +251,7 @@ hm_w = heatmap!(ax_w, xw, zw, wₙ; colorrange = (-wlim, wlim), colormap = :bala Colorbar(fig[3, 2], hm_w, label = "m s⁻¹") ax_N² = Axis(fig[4, 1]; title = "stratification N²", axis_kwargs...) -hm_N² = heatmap!(ax_N², xN², zN², N²ₙ; colorrange = (0.95Nᵢ², 1.05Nᵢ²), colormap = :thermal) +hm_N² = heatmap!(ax_N², xN², zN², N²ₙ; colorrange = (0.9Nᵢ², 1.1Nᵢ²), colormap = :thermal) Colorbar(fig[4, 2], hm_N², label = "s⁻²") fig From 1b3950749af5b1ff8f4709e93c560f74f8cf89a4 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Mon, 5 Jun 2023 22:03:31 +1000 Subject: [PATCH 22/43] define H --- examples/internal_tide.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 6b6b71a752..b0617a2197 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -22,9 +22,11 @@ using Oceananigans.Units Nx, Nz = 250, 125 +H = 2kilometers + underlying_grid = RectilinearGrid(size = (Nx, Nz), x = (-1000kilometers, 1000kilometers), - z = (-2kilometers, 0), + z = (-H, 0), halo = (4, 4), topology = (Periodic, Flat, Bounded)) From 89fd8536fbed6caeae10ce42dcb2666ec7135537 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 Jun 2023 06:34:56 +1000 Subject: [PATCH 23/43] undo pcg modifications --- .../pcg_implicit_free_surface_solver.jl | 35 +++++++++---------- 1 file changed, 16 insertions(+), 19 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl b/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl index 9bb5913524..cfa612d718 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/pcg_implicit_free_surface_solver.jl @@ -2,7 +2,7 @@ using Oceananigans.Solvers using Oceananigans.Operators using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBottom using Oceananigans.Architectures -using Oceananigans.Grids: with_halo, isrectilinear, halo_size +using Oceananigans.Grids: with_halo, isrectilinear using Oceananigans.Architectures: device import Oceananigans.Solvers: solve!, precondition! @@ -43,8 +43,8 @@ step `Δt`, gravitational acceleration `g`, and free surface at time-step `n` ` """ function PCGImplicitFreeSurfaceSolver(grid::AbstractGrid, settings, gravitational_acceleration=nothing) # Initialize vertically integrated lateral face areas - ∫ᶻ_Axᶠᶜᶜ = Field((Face, Center, Nothing), grid) - ∫ᶻ_Ayᶜᶠᶜ = Field((Center, Face, Nothing), grid) + ∫ᶻ_Axᶠᶜᶜ = Field((Face, Center, Nothing), with_halo((3, 3, 1), grid)) + ∫ᶻ_Ayᶜᶠᶜ = Field((Center, Face, Nothing), with_halo((3, 3, 1), grid)) vertically_integrated_lateral_areas = (xᶠᶜᶜ = ∫ᶻ_Axᶠᶜᶜ, yᶜᶠᶜ = ∫ᶻ_Ayᶜᶠᶜ) @@ -82,10 +82,10 @@ function solve!(η, implicit_free_surface_solver::PCGImplicitFreeSurfaceSolver, # Take explicit step first? We haven't found improvement from this yet, but perhaps it will # help eventually. #explicit_ab2_step_free_surface!(free_surface, model, Δt, χ) - + ∫ᶻA = implicit_free_surface_solver.vertically_integrated_lateral_areas solver = implicit_free_surface_solver.preconditioned_conjugate_gradient_solver - + # solve!(x, solver, b, args...) solves A*x = b for x. solve!(η, solver, rhs, ∫ᶻA.xᶠᶜᶜ, ∫ᶻA.yᶜᶠᶜ, g, Δt) @@ -118,9 +118,7 @@ end end """ - implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) - -Return `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n` +Returns `L(ηⁿ)`, where `ηⁿ` is the free surface displacement at time step `n` and `L` is the linear operator that arises in an implicit time step for the free surface displacement `η`. @@ -144,7 +142,7 @@ end """ _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) -Return the left side of the "implicit ``η`` equation" +Return the left side of the "implicit η equation" ```math (∇ʰ⋅ H ∇ʰ - 1 / (g Δt²)) ηⁿ⁺¹ = 1 / (g Δt) ∇ʰ ⋅ Q★ - 1 / (g Δt²) ηⁿ @@ -159,7 +157,7 @@ and arranged in a symmetric form by multiplying by horizontal areas Az: δⁱÂʷ∂ˣηⁿ⁺¹ + δʲÂˢ∂ʸηⁿ⁺¹ - Az ηⁿ⁺¹ / (g Δt²) = 1 / (g Δt) (δⁱÂʷu̅ˢᵗᵃʳ + δʲÂˢv̅ˢᵗᵃʳ) - Az ηⁿ / (g Δt²) ``` -where ̂ indicates a vertical integral, and +where ̂ indicates a vertical integral, and ̅ indicates a vertical average """ @kernel function _implicit_free_surface_linear_operation!(L_ηⁿ⁺¹, grid, ηⁿ⁺¹, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) @@ -174,19 +172,20 @@ end ##### """ -Add `- H⁻¹ ∇H ⋅ ∇ηⁿ` to the right-hand-side. +add to the rhs - H⁻¹ ∇H ⋅ ∇ηⁿ to the rhs... """ @inline function precondition!(P_r, preconditioner::FFTImplicitFreeSurfaceSolver, r, η, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) poisson_solver = preconditioner.fft_poisson_solver arch = architecture(poisson_solver) grid = preconditioner.three_dimensional_grid Az = grid.Δxᶜᵃᵃ * grid.Δyᵃᶜᵃ # assume horizontal regularity - Lz = grid.Lz + Lz = grid.Lz launch!(arch, grid, :xy, fft_preconditioner_right_hand_side!, poisson_solver.storage, r, η, grid, Az, Lz) + return solve!(P_r, preconditioner, poisson_solver.storage, g, Δt) end @@ -232,22 +231,20 @@ struct DiagonallyDominantInversePreconditioner end """ _diagonally_dominant_precondition!(P_r, grid, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) -Return the diagonally dominant inverse preconditioner applied to the residuals consistently -with `M = D⁻¹(I - (A - D)D⁻¹) ≈ A⁻¹` where `I` is the identity matrix, `A` is the linear -operator applied to the free surface `η`, and `D` is the diagonal of `A`. +Return the diagonally dominant inverse preconditioner applied to the residuals consistently with + `M = D⁻¹(I - (A - D)D⁻¹) ≈ A⁻¹` where `I` is the Identity matrix, +A is the linear operator applied to η, and D is the diagonal of A. ```math P_r = M * r ``` - which expanded in components is - ```math P_rᵢⱼ = rᵢⱼ / Acᵢⱼ - 1 / Acᵢⱼ ( Ax⁻ / Acᵢ₋₁ rᵢ₋₁ⱼ + Ax⁺ / Acᵢ₊₁ rᵢ₊₁ⱼ + Ay⁻ / Acⱼ₋₁ rᵢⱼ₋₁+ Ay⁺ / Acⱼ₊₁ rᵢⱼ₊₁ ) ``` -where `Ac`, `Ax⁻`, `Ax⁺`, `Ay⁻` and `Ay⁺` are the coefficients of `ηᵢⱼ`, `ηᵢ₋₁ⱼ`, `ηᵢ₊₁ⱼ`, `ηᵢⱼ₋₁`, -and `ηᵢⱼ₊₁` in `_implicit_free_surface_linear_operation!` +where `Ac`, `Ax⁻`, `Ax⁺`, `Ay⁻` and `Ay⁺` are the coefficients of +`ηᵢⱼ`, `ηᵢ₋₁ⱼ`, `ηᵢ₊₁ⱼ`, `ηᵢⱼ₋₁` and `ηᵢⱼ₊₁` in `_implicit_free_surface_linear_operation!` """ function diagonally_dominant_precondition!(P_r, r, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt) grid = ∫ᶻ_Axᶠᶜᶜ.grid From b2eb120fc9017c7c7f9b4c6383990428181eae4d Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 Jun 2023 06:35:48 +1000 Subject: [PATCH 24/43] bring back all examples --- docs/make.jl | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 3b2d62160c..b33ef026a6 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -28,18 +28,18 @@ const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples") const OUTPUT_DIR = joinpath(@__DIR__, "src/generated") examples = [ - # "One-dimensional diffusion" => "one_dimensional_diffusion", - # "Two-dimensional turbulence" => "two_dimensional_turbulence", - # "Internal wave" => "internal_wave", + "One-dimensional diffusion" => "one_dimensional_diffusion", + "Two-dimensional turbulence" => "two_dimensional_turbulence", + "Internal wave" => "internal_wave", "Internal tide by a seamount" => "internal_tide", - # "Convecting plankton" => "convecting_plankton", - # "Ocean wind mixing and convection" => "ocean_wind_mixing_and_convection", - # "Langmuir turbulence" => "langmuir_turbulence", - # "Baroclinic adjustment" => "baroclinic_adjustment", - # "Kelvin-Helmholtz instability" => "kelvin_helmholtz_instability", - # "Shallow water Bickley jet" => "shallow_water_Bickley_jet", - # "Horizontal convection" => "horizontal_convection", - # "Tilted bottom boundary layer" => "tilted_bottom_boundary_layer" + "Convecting plankton" => "convecting_plankton", + "Ocean wind mixing and convection" => "ocean_wind_mixing_and_convection", + "Langmuir turbulence" => "langmuir_turbulence", + "Baroclinic adjustment" => "baroclinic_adjustment", + "Kelvin-Helmholtz instability" => "kelvin_helmholtz_instability", + "Shallow water Bickley jet" => "shallow_water_Bickley_jet", + "Horizontal convection" => "horizontal_convection", + "Tilted bottom boundary layer" => "tilted_bottom_boundary_layer" ] example_scripts = [ filename * ".jl" for (title, filename) in examples ] @@ -144,10 +144,10 @@ makedocs(bib, format = format, pages = pages, modules = [Oceananigans], - doctest = false, + doctest = true, strict = false, clean = true, - checkdocs = :none + checkdocs = :exports ) @info "Clean up temporary .jld2 and .nc output created by doctests or literated examples..." From 273fd1e00e2209fc38d1b3029d84760816cbf0e5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 Jun 2023 06:36:27 +1000 Subject: [PATCH 25/43] bring back strictness --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index b33ef026a6..c1ff86c9c9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -145,7 +145,7 @@ makedocs(bib, pages = pages, modules = [Oceananigans], doctest = true, - strict = false, + strict = true, clean = true, checkdocs = :exports ) From 7810d49a08db55713f482fa143c9721b00d7c83e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 Jun 2023 06:37:29 +1000 Subject: [PATCH 26/43] simpler forcing description --- examples/internal_tide.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index b0617a2197..69ae3b54ec 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -65,10 +65,10 @@ lines!(ax, xC/1e3, bottom_boundary) fig -# Now we want to add a barotropic tide forcing. For example, to add the ``M_2`` tidal forcing +# Now we want to add a barotropic tide forcing. For example, to add the lunar semi-diurnal ``M_2`` tide # we need to add forcing in the ``u``-momentum equation of the form: # ```math -# \partial_t u = \dotsb + F_0 \sin(\omega_2 t) +# F_0 \sin(\omega_2 t) # ``` # where ``\omega_2 = 2π / T_2``, with ``T_2 = 12.421 \,\mathrm{hours}`` the period of the ``M_2`` tide. From 8aa442da698e398183e0d2d8c2138bb5e2dfcded Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 6 Jun 2023 06:39:03 +1000 Subject: [PATCH 27/43] \varepsilon -> \epsilon --- examples/internal_tide.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 69ae3b54ec..12228cdf6b 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -101,9 +101,9 @@ coriolis = FPlane(latitude = -45) T₂ = 12.421hours ω₂ = 2π / T₂ # radians/sec -ε = 0.1 # excursion parameter +ϵ = 0.1 # excursion parameter -U_tidal = ε * ω₂ * width +U_tidal = ϵ * ω₂ * width tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂ From f13a33ef8e43a888e00b7063a231e117c5b482bd Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 7 Jun 2023 08:50:26 +1000 Subject: [PATCH 28/43] code alignment --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index c1ff86c9c9..ac3752b25d 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -31,7 +31,7 @@ examples = [ "One-dimensional diffusion" => "one_dimensional_diffusion", "Two-dimensional turbulence" => "two_dimensional_turbulence", "Internal wave" => "internal_wave", - "Internal tide by a seamount" => "internal_tide", + "Internal tide by a seamount" => "internal_tide", "Convecting plankton" => "convecting_plankton", "Ocean wind mixing and convection" => "ocean_wind_mixing_and_convection", "Langmuir turbulence" => "langmuir_turbulence", From 0ba8f27b9555aaaa82666fb5a013bed469a210d5 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 30 Jun 2023 11:28:51 +1000 Subject: [PATCH 29/43] order according to most computationally expensive --- docs/make.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 83510aa5a1..436508abb3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -46,15 +46,16 @@ end "baroclinic_adjustment.jl", "kelvin_helmholtz_instability.jl", "langmuir_turbulence.jl", - "ocean_wind_mixing_and_convection.jl", "horizontal_convection.jl", - "convecting_plankton.jl", "tilted_bottom_boundary_layer.jl", + "ocean_wind_mixing_and_convection.jl", + "convecting_plankton.jl", "shallow_water_Bickley_jet.jl", "internal_tide", "two_dimensional_turbulence.jl", "internal_wave.jl", "one_dimensional_diffusion.jl", + ] @info string("Executing the examples using ", Distributed.nprocs(), " processes") end From 4f051c6a1cb9c0dd2cc1cac1868c5f49266ca7d9 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 30 Jun 2023 14:27:56 +1000 Subject: [PATCH 30/43] swe example takes long time --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 436508abb3..6d2700fd61 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -43,6 +43,7 @@ end # The examples that take longer to run should be first. This ensures thats # docs built using extra workers is as efficient as possible. example_scripts = [ + "shallow_water_Bickley_jet.jl", "baroclinic_adjustment.jl", "kelvin_helmholtz_instability.jl", "langmuir_turbulence.jl", @@ -50,7 +51,6 @@ end "tilted_bottom_boundary_layer.jl", "ocean_wind_mixing_and_convection.jl", "convecting_plankton.jl", - "shallow_water_Bickley_jet.jl", "internal_tide", "two_dimensional_turbulence.jl", "internal_wave.jl", From 71d2661c0259577b38b6cd412981249b8e9e3b9f Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 16 Jan 2024 13:36:16 +0200 Subject: [PATCH 31/43] updates deprecated syntax --- examples/internal_tide.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 12228cdf6b..50d3e53d60 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -35,16 +35,16 @@ underlying_grid = RectilinearGrid(size = (Nx, Nz), # it returns the ``z`` of the bottom, or *(ii)* a two-dimensional array with the values of ``z`` at # the bottom cell centers. # -# In this example we'd like to have a small Gaussian hill at the center of the domain. +# In this example we'd like to have a Gaussian hill at the center of the domain. # # ```math # h(x) = -H + h_0 \exp(-x^2 / 2σ^2) # ``` -h₀ = 250 # m +h₀ = 250meters width = 20kilometers hill(x) = h₀ * exp(-x^2 / 2width^2) -bottom(x, y) = - H + hill(x) +bottom(x) = - H + hill(x) grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom)) @@ -59,7 +59,7 @@ fig = Figure(resolution = (700, 200)) ax = Axis(fig[1, 1], xlabel="x [km]", ylabel="z [m]", - limits=((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz, -4grid.Lz/5))) + limits=((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz, 0))) lines!(ax, xC/1e3, bottom_boundary) From 4b7342815562f862006043753adcd033a55e6eee Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 16 Jan 2024 14:27:46 +0200 Subject: [PATCH 32/43] minor code organization --- src/OutputWriters/fetch_output.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OutputWriters/fetch_output.jl b/src/OutputWriters/fetch_output.jl index 768f8d3514..5553fb2dd8 100644 --- a/src/OutputWriters/fetch_output.jl +++ b/src/OutputWriters/fetch_output.jl @@ -1,6 +1,7 @@ using CUDA using Oceananigans.Fields: AbstractField, compute_at!, ZeroField +using Oceananigans.ImmersedBoundaries: mask_immersed_field! using Oceananigans.Models.LagrangianParticleTracking: LagrangianParticles # Needed to support `fetch_output` with `model::Nothing`. @@ -43,4 +44,3 @@ function fetch_and_convert_output(output, model, writer) end fetch_and_convert_output(output::ZeroField, model, writer) = zero(eltype(output)) - From 0c2201244ecd4a87d75a1a053eec5c4bc3cdd25a Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 16 Jan 2024 14:28:03 +0200 Subject: [PATCH 33/43] remove stray line --- src/ImmersedBoundaries/mask_immersed_field.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/ImmersedBoundaries/mask_immersed_field.jl b/src/ImmersedBoundaries/mask_immersed_field.jl index 4bd2878eec..8a04082d89 100644 --- a/src/ImmersedBoundaries/mask_immersed_field.jl +++ b/src/ImmersedBoundaries/mask_immersed_field.jl @@ -56,4 +56,3 @@ end @inline scalar_mask(i, j, k, grid, ::AbstractGridFittedBoundary, LX, LY, LZ, value, field, mask=peripheral_node) = @inbounds ifelse(mask(i, j, k, grid, LX, LY, LZ), value, field[i, j, k]) - From 398388286b9a51579cbc29a13ba65b0ca52c11df Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Tue, 16 Jan 2024 14:28:37 +0200 Subject: [PATCH 34/43] update deprecated syntax --- examples/internal_tide.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 50d3e53d60..394521ed25 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -97,7 +97,7 @@ coriolis = FPlane(latitude = -45) # Now we have everything we require to construct the tidal forcing given a value of the # excursion parameter. -# + T₂ = 12.421hours ω₂ = 2π / T₂ # radians/sec @@ -107,7 +107,7 @@ U_tidal = ϵ * ω₂ * width tidal_forcing_amplitude = U_tidal * (ω₂^2 - coriolis.f^2) / ω₂ -@inline tidal_forcing(x, y, z, t, p) = p.tidal_forcing_amplitude * sin(p.ω₂ * t) +@inline tidal_forcing(x, z, t, p) = p.tidal_forcing_amplitude * sin(p.ω₂ * t) u_forcing = Forcing(tidal_forcing, parameters=(; tidal_forcing_amplitude, ω₂)) @@ -116,8 +116,7 @@ u_forcing = Forcing(tidal_forcing, parameters=(; tidal_forcing_amplitude, ω₂) # We built a `HydrostaticFreeSurfaceModel` with a `SplitExplicitFreeSurface` solver. # We limit our maximum timestep to 3 minutes. -max_Δt = 5minutes -free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7, max_Δt) +free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7) model = HydrostaticFreeSurfaceModel(; grid, coriolis, free_surface, buoyancy = BuoyancyTracer(), @@ -128,16 +127,16 @@ model = HydrostaticFreeSurfaceModel(; grid, coriolis, free_surface, # We initialize the model with the tidal flow and a linear stratification. -uᵢ(x, y, z) = U_tidal +uᵢ(x, z) = U_tidal Nᵢ² = 1e-4 # [s⁻²] initial buoyancy frequency / stratification -bᵢ(x, y, z) = Nᵢ² * z +bᵢ(x, z) = Nᵢ² * z set!(model, u=uᵢ, b=bᵢ) # Now let's built a `Simulation`. -Δt = max_Δt +Δt = 5minutes stop_time = 4days simulation = Simulation(model; Δt, stop_time) From 135056919f1fcf5365dba61c114d898e0005b0fe Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 13:58:53 +0200 Subject: [PATCH 35/43] update internal_tide example --- examples/internal_tide.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 394521ed25..76ebb8a7f5 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -48,20 +48,20 @@ bottom(x) = - H + hill(x) grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom)) -# The hill is small; here's how it looks (note that we don't plot all the way to the ocean surface). +# Let's see how the domain with the bathymetry looks. -xC = xnodes(grid, Center()) -bottom_boundary = grid.immersed_boundary.bottom_height[1:Nx, 1] +x = xnodes(grid, Center()) +bottom_boundary = interior(grid.immersed_boundary.bottom_height, :, 1, 1) using CairoMakie -fig = Figure(resolution = (700, 200)) +fig = Figure(size = (700, 200)) ax = Axis(fig[1, 1], xlabel="x [km]", ylabel="z [m]", limits=((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz, 0))) -lines!(ax, xC/1e3, bottom_boundary) +lines!(ax, x/1e3, bottom_boundary) fig @@ -147,8 +147,7 @@ using Printf wall_clock = Ref(time_ns()) -function print_progress(sim) - +function progress(sim) elapsed = 1e-9 * (time_ns() - wall_clock[]) msg = @sprintf("iteration: %d, time: %s, wall time: %s, max|w|: %6.3e, m s⁻¹\n", @@ -162,7 +161,7 @@ function print_progress(sim) return nothing end -simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(200)) +add_callback!(simulation, progress, name=:progress, IterationInterval(200)) nothing #hide # ## Diagnostics/Output @@ -239,7 +238,7 @@ axis_kwargs = (xlabel = "x [km]", ulim = maximum(abs, u′_t[end]) wlim = maximum(abs, w_t[end]) -fig = Figure(resolution = (700, 900)) +fig = Figure(size = (700, 900)) fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) From 979c2a7cd10c1f2054cad996622cafd2e3fb3ffa Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 14:26:25 +0200 Subject: [PATCH 36/43] mask immersed output with 0 --- src/OutputWriters/fetch_output.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/OutputWriters/fetch_output.jl b/src/OutputWriters/fetch_output.jl index 43c333dfe7..1b02cc6849 100644 --- a/src/OutputWriters/fetch_output.jl +++ b/src/OutputWriters/fetch_output.jl @@ -17,6 +17,7 @@ fetch_output(output, model) = output(model) function fetch_output(field::AbstractField, model) compute_at!(field, time(model)) + mask_immersed_field!(field) return parent(field) end From 06e20a54941fa073d739f2fe5eb155405bffe66e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 18:00:57 +0200 Subject: [PATCH 37/43] revert change --- src/OutputWriters/fetch_output.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/OutputWriters/fetch_output.jl b/src/OutputWriters/fetch_output.jl index 1b02cc6849..a86d1d6f9b 100644 --- a/src/OutputWriters/fetch_output.jl +++ b/src/OutputWriters/fetch_output.jl @@ -1,7 +1,6 @@ using CUDA using Oceananigans.Fields: AbstractField, compute_at!, ZeroField -using Oceananigans.ImmersedBoundaries: mask_immersed_field! using Oceananigans.Models.LagrangianParticleTracking: LagrangianParticles # TODO: figure out how to support this @@ -17,7 +16,6 @@ fetch_output(output, model) = output(model) function fetch_output(field::AbstractField, model) compute_at!(field, time(model)) - mask_immersed_field!(field) return parent(field) end From da4502a0843ce266be6f755868d30be5a6913e50 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 18:01:18 +0200 Subject: [PATCH 38/43] add post-process masking step --- examples/internal_tide.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 76ebb8a7f5..921ceff525 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -204,6 +204,16 @@ N²_t = FieldTimeSeries(saved_output_filename, "N²") times = u′_t.times nothing #hide +# For visualization purposes, we mask the region below the bathymetry with NaNs. + +using Oceananigans.ImmersedBoundaries: mask_immersed_field! + +for n in 1:length(times) + mask_immersed_field!(u′_t[n], NaN) + mask_immersed_field!(w_t[n], NaN) + mask_immersed_field!(N²_t[n], NaN) +end + # We retrieve each field's coordinates and convert from meters to kilometers. xu, yu, zu = nodes(u′_t[1]) ./ 1e3 From 84bebe8e2407aa3f6a385f42509a22c71f84b7d3 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 18:02:18 +0200 Subject: [PATCH 39/43] revert changes --- src/ImmersedBoundaries/mask_immersed_field.jl | 1 + src/OutputWriters/fetch_output.jl | 1 + 2 files changed, 2 insertions(+) diff --git a/src/ImmersedBoundaries/mask_immersed_field.jl b/src/ImmersedBoundaries/mask_immersed_field.jl index 8a04082d89..4bd2878eec 100644 --- a/src/ImmersedBoundaries/mask_immersed_field.jl +++ b/src/ImmersedBoundaries/mask_immersed_field.jl @@ -56,3 +56,4 @@ end @inline scalar_mask(i, j, k, grid, ::AbstractGridFittedBoundary, LX, LY, LZ, value, field, mask=peripheral_node) = @inbounds ifelse(mask(i, j, k, grid, LX, LY, LZ), value, field[i, j, k]) + diff --git a/src/OutputWriters/fetch_output.jl b/src/OutputWriters/fetch_output.jl index a86d1d6f9b..99de6cc505 100644 --- a/src/OutputWriters/fetch_output.jl +++ b/src/OutputWriters/fetch_output.jl @@ -48,3 +48,4 @@ function fetch_and_convert_output(output, model, writer) end fetch_and_convert_output(output::ZeroField, model, writer) = zero(eltype(output)) + From cf18546303c3c619ac66630a32525a8a474a91c7 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 18:04:35 +0200 Subject: [PATCH 40/43] find max abs values before masking --- examples/internal_tide.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 921ceff525..d6f42d2a94 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -201,6 +201,9 @@ u′_t = FieldTimeSeries(saved_output_filename, "u′") w_t = FieldTimeSeries(saved_output_filename, "w") N²_t = FieldTimeSeries(saved_output_filename, "N²") +umax = maximum(abs, u′_t[end]) +wmax = maximum(abs, w_t[end]) + times = u′_t.times nothing #hide @@ -245,19 +248,16 @@ axis_kwargs = (xlabel = "x [km]", limits = ((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz/1e3, 0)), # note conversion to kilometers titlesize = 20) -ulim = maximum(abs, u′_t[end]) -wlim = maximum(abs, w_t[end]) - fig = Figure(size = (700, 900)) fig[1, :] = Label(fig, title, fontsize=24, tellwidth=false) ax_u = Axis(fig[2, 1]; title = "u'-velocity", axis_kwargs...) -hm_u = heatmap!(ax_u, xu, zu, u′ₙ; colorrange = (-ulim, ulim), colormap = :balance) +hm_u = heatmap!(ax_u, xu, zu, u′ₙ; colorrange = (-umax, umax), colormap = :balance) Colorbar(fig[2, 2], hm_u, label = "m s⁻¹") ax_w = Axis(fig[3, 1]; title = "w-velocity", axis_kwargs...) -hm_w = heatmap!(ax_w, xw, zw, wₙ; colorrange = (-wlim, wlim), colormap = :balance) +hm_w = heatmap!(ax_w, xw, zw, wₙ; colorrange = (-wmax, wmax), colormap = :balance) Colorbar(fig[3, 2], hm_w, label = "m s⁻¹") ax_N² = Axis(fig[4, 1]; title = "stratification N²", axis_kwargs...) From cf1ad7b6f3231e5c633e5f0cdc728ce85876db52 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Fri, 8 Mar 2024 10:56:17 +0200 Subject: [PATCH 41/43] plot domain with band --- examples/internal_tide.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index d6f42d2a94..c6b74a7eb5 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -48,10 +48,11 @@ bottom(x) = - H + hill(x) grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom)) -# Let's see how the domain with the bathymetry looks. +# Let's see how the domain with the bathymetry is. x = xnodes(grid, Center()) bottom_boundary = interior(grid.immersed_boundary.bottom_height, :, 1, 1) +top_boundary = 0*x using CairoMakie @@ -61,7 +62,7 @@ ax = Axis(fig[1, 1], ylabel="z [m]", limits=((-grid.Lx/2e3, grid.Lx/2e3), (-grid.Lz, 0))) -lines!(ax, x/1e3, bottom_boundary) +band!(ax, x/1e3, bottom_boundary, top_boundary, color = :mediumblue) fig @@ -86,7 +87,7 @@ fig # the tidal velocity and tidal forcing amplitudes are related via: # # ```math -# U_{\mathrm{tidal}} = \frac{\omega_2}{\omega_2^2 - f^2}F_0 +# U_{\mathrm{tidal}} = \frac{\omega_2}{\omega_2^2 - f^2} F_0 # ``` # # Now we have the way to find the value of the tidal forcing amplitude that would correspond to a From ca437145013e82380bbbc8b1d700656eeba78607 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sat, 9 Mar 2024 14:23:36 +0200 Subject: [PATCH 42/43] simpler masking --- examples/internal_tide.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index c6b74a7eb5..8209e80d70 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -212,10 +212,8 @@ nothing #hide using Oceananigans.ImmersedBoundaries: mask_immersed_field! -for n in 1:length(times) - mask_immersed_field!(u′_t[n], NaN) - mask_immersed_field!(w_t[n], NaN) - mask_immersed_field!(N²_t[n], NaN) +for φ_t in (u′_t, w_t, N²_t), n in 1:length(times) + mask_immersed_field!(φ_t[n], NaN) end # We retrieve each field's coordinates and convert from meters to kilometers. From 6777dd7e01d44b2b7fdb70e0481898aa6be7c306 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Sun, 10 Mar 2024 21:51:12 +0200 Subject: [PATCH 43/43] fix missing .jl --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 5c10255575..568ffe94ec 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -49,7 +49,7 @@ end "langmuir_turbulence.jl", "ocean_wind_mixing_and_convection.jl", "horizontal_convection.jl", - "internal_tide", + "internal_tide.jl", "convecting_plankton.jl", "tilted_bottom_boundary_layer.jl", "two_dimensional_turbulence.jl",