diff --git a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl index a5bc6435c7..b91ac36249 100644 --- a/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl +++ b/src/BoundaryConditions/flat_extrapolation_open_boundary_matching_scheme.jl @@ -128,7 +128,7 @@ end spacing_factor = Δz₁ / (Δz₂ + Δz₃) - gradient_free_ϕ = @inbounds ϕ[i, j, 3] - (ϕ[i, k, 2] - ϕ[i, j, 4]) * spacing_factor + gradient_free_ϕ = @inbounds ϕ[i, j, 3] - (ϕ[i, j, 2] - ϕ[i, j, 4]) * spacing_factor @inbounds ϕ[i, j, 1] = relax(i, j, grid, gradient_free_ϕ, bc, clock, model_fields) @@ -149,4 +149,4 @@ end @inbounds ϕ[i, j, k] = relax(i, j, grid, gradient_free_ϕ, bc, clock, model_fields) return nothing -end \ No newline at end of file +end diff --git a/src/Grids/orthogonal_spherical_shell_grid.jl b/src/Grids/orthogonal_spherical_shell_grid.jl index 1567b2a5fb..7d6edd4f5b 100644 --- a/src/Grids/orthogonal_spherical_shell_grid.jl +++ b/src/Grids/orthogonal_spherical_shell_grid.jl @@ -76,7 +76,7 @@ end const OSSG = OrthogonalSphericalShellGrid const ZRegOSSG = OrthogonalSphericalShellGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Number} const ZRegOrthogonalSphericalShellGrid = ZRegOSSG -const ConformalCubedSpherePanel = OrthogonalSphericalShellGrid{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:CubedSphereConformalMapping} +const ConformalCubedSpherePanel = OrthogonalSphericalShellGrid{<:Any, FullyConnected, FullyConnected, <:Any, <:Any, <:Any, <:Any, <:CubedSphereConformalMapping} # convenience constructor for OSSG without any conformal_mapping properties OrthogonalSphericalShellGrid(architecture, Nx, Ny, Nz, Hx, Hy, Hz, Lz, diff --git a/src/Operators/vorticity_operators.jl b/src/Operators/vorticity_operators.jl index 170a013a28..bf1f95908d 100644 --- a/src/Operators/vorticity_operators.jl +++ b/src/Operators/vorticity_operators.jl @@ -10,6 +10,12 @@ The vertical vorticity associated with horizontal velocities ``u`` and ``v``. """ @inline ζ₃ᶠᶠᶜ(i, j, k, grid, u, v) = Γᶠᶠᶜ(i, j, k, grid, u, v) / Azᶠᶠᶜ(i, j, k, grid) +# South-west, south-east, north-west, north-east corners +@inline on_south_west_corner(i, j, grid) = (i == 1) & (j == 1) +@inline on_south_east_corner(i, j, grid) = (i == grid.Nx+1) & (j == 1) +@inline on_north_east_corner(i, j, grid) = (i == grid.Nx+1) & (j == grid.Ny+1) +@inline on_north_west_corner(i, j, grid) = (i == 1) & (j == grid.Ny+1) + ##### ##### Vertical circulation at the corners of the cubed sphere needs to treated in a special manner. ##### See: https://github.com/CliMA/Oceananigans.jl/issues/1584 @@ -21,24 +27,13 @@ The vertical vorticity associated with horizontal velocities ``u`` and ``v``. The vertical circulation associated with horizontal velocities ``u`` and ``v``. """ @inline function Γᶠᶠᶜ(i, j, k, grid::ConformalCubedSpherePanel, u, v) - # South-west corner - if i == 1 && j == 1 - return Δy_qᶜᶠᶜ(i, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u) - - # South-east corner - elseif i == grid.Nx+1 && j == 1 - return - Δy_qᶜᶠᶜ(i-1, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u) - - # North-west corner - elseif i == 1 && j == grid.Ny+1 - return Δy_qᶜᶠᶜ(i, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u) - - # North-east corner - elseif i == grid.Nx+1 && j == grid.Ny+1 - return - Δy_qᶜᶠᶜ(i-1, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u) - - # Not a corner - else - return δxᶠᵃᵃ(i, j, k, grid, Δy_qᶜᶠᶜ, v) - δyᵃᶠᵃ(i, j, k, grid, Δx_qᶠᶜᶜ, u) - end + Hx, Hy = grid.Hx, grid.Hy + Γ = ifelse(on_south_west_corner(i, j, grid) | on_north_west_corner(i, j, grid), + Δy_qᶜᶠᶜ(i, j, k, grid, v) - Δx_qᶠᶜᶜ(i, j, k, grid, u) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u), + ifelse(on_south_east_corner(i, j, grid) | on_north_east_corner(i, j, grid), + - Δy_qᶜᶠᶜ(i-1, j, k, grid, v) + Δx_qᶠᶜᶜ(i, j-1, k, grid, u) - Δx_qᶠᶜᶜ(i, j, k, grid, u), + δxᶠᶠᶜ(i, j, k, grid, Δy_qᶜᶠᶜ, v) - δyᶠᶠᶜ(i, j, k, grid, Δx_qᶠᶜᶜ, u) + ) + ) + return Γ end diff --git a/validation/open_boundaries/oscillating_flow.jl b/validation/open_boundaries/oscillating_flow.jl index 6bf48337ea..a2e394d6c1 100644 --- a/validation/open_boundaries/oscillating_flow.jl +++ b/validation/open_boundaries/oscillating_flow.jl @@ -1,12 +1,11 @@ -# This validation script shows open boundaries working in a simple case where the -# oscillates sinusoidally so changes sign across two open boundaries. This is similar -# to a more realistic case where we know some arbitary external conditions. -# This necessitates using a combination allowing information to exit the domain, in -# this case by setting the wall normal velocity gradient to zero, as well as forcing -# to the external value in this example by relaxing to it. - -# This case also has a stretched grid to validate the zero wall normal velocity -# gradient matching scheme on a stretched grid. +# This validation script shows open boundaries working in a simple case where a flow past a 2D +# cylinder oscillates in two directions. All boundaries have the same +# `FlatExtrapolationOpenBoundaryCondition`s. This is similar to a more realistic case where we know +# some arbitary external conditions. First we test an xy flow and then we test an xz flow (the +# forcings and boundary conditions originally designed for `v` aere then used for `w` without +# modification). +# +# This case also has a stretched grid to validate the matching scheme on a stretched grid. using Oceananigans, CairoMakie using Oceananigans.BoundaryConditions: FlatExtrapolationOpenBoundaryCondition @@ -22,133 +21,125 @@ end architecture = CPU() # model parameters -Re = 200 U = 1 -D = 1. -resolution = D / 10 - -# add extra downstream distance to see if the solution near the cylinder changes -extra_downstream = 0 +D = 1.0 +T = 50 cylinder = Cylinder(; D) -x = (-5, 5 + extra_downstream) .* D - -Ny = Int(10 / resolution) - -Nx = Int((10 + extra_downstream) / resolution) - -function Δy(j) - if Ny/2 - 2/resolution < j < Ny/2 + 2/resolution - return resolution - elseif j <= Ny/2 - 2/resolution - return resolution * (1 + (Ny/2 - 2/resolution - j) / (Ny/2 - 2/resolution)) - elseif j >= Ny/2 + 2/resolution - return resolution * (1 + (j - Ny/2 - 2/resolution) / (Ny/2 - 2/resolution)) - else - Throw(ArgumentError("$j not in range")) - end -end - -y(j) = sum(Δy.([1:j;])) - sum(Δy.([1:Ny;]))/2 - -ν = U * D / Re - -closure = ScalarDiffusivity(;ν, κ = ν) - -grid = RectilinearGrid(architecture; topology = (Bounded, Bounded, Flat), size = (Nx, Ny), x = y, y = x) - -T = 20 / U - -@inline u(t, p) = p.U * sin(t * 2π / p.T) -@inline u(y, t, p) = u(t, p) - -relaxation_timescale = 0.15 - -u_boundaries = FieldBoundaryConditions(east = FlatExtrapolationOpenBoundaryCondition(u; relaxation_timescale, parameters = (; U, T)), - west = FlatExtrapolationOpenBoundaryCondition(u; relaxation_timescale, parameters = (; U, T)), - south = GradientBoundaryCondition(0), - north = GradientBoundaryCondition(0)) - -v_boundaries = FieldBoundaryConditions(east = GradientBoundaryCondition(0), - west = GradientBoundaryCondition(0), - south = FlatExtrapolationOpenBoundaryCondition(0; relaxation_timescale), - north = FlatExtrapolationOpenBoundaryCondition(0; relaxation_timescale)) +L = 10 +Nx = Ny = 40 -Δt = .3 * resolution / U +β = 0.2 +x_faces(i) = L/2 * (β * ((2 * (i - 1)) / Nx - 1)^3 + (2 * (i - 1)) / Nx - 1) / (β + 1) +y = (-L/2, +L/2) .* D -@show Δt +xygrid = RectilinearGrid(architecture; topology = (Bounded, Bounded, Flat), size = (Nx, Ny), x = x_faces, y) +xzgrid = RectilinearGrid(architecture; topology = (Bounded, Flat, Bounded), size = (Nx, Ny), x = x_faces, z = y) -u_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder) -v_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder) +Δt = .5 * minimum_xspacing(xygrid) / abs(U) -model = NonhydrostaticModel(; grid, - closure, - forcing = (u = u_forcing, v = v_forcing), - boundary_conditions = (u = u_boundaries, v = v_boundaries)) +@inline u∞(y, t, p) = p.U * cos(t * 2π / p.T) * (1 + 0.01 * randn()) +@inline v∞(x, t, p) = p.U * sin(t * 2π / p.T) * (1 + 0.01 * randn()) -@info "Constructed model" +function run_cylinder(grid, boundary_conditions; plot=true, stop_time = 50, simname = "") + @info "Testing $simname with grid" grid -# initial noise to induce turbulance faster -set!(model, u = (x, y) -> randn() * U * 0.01, v = (x, y) -> randn() * U * 0.01) + cylinder_forcing = Relaxation(; rate = 1 / (2 * Δt), mask = cylinder) -@info "Set initial conditions" + global model = NonhydrostaticModel(; grid, + advection = UpwindBiasedFifthOrder(), + forcing = (u = cylinder_forcing, v = cylinder_forcing, w = cylinder_forcing), + boundary_conditions) -simulation = Simulation(model; Δt = Δt, stop_time = 300) + @info "Constructed model" -wizard = TimeStepWizard(cfl = 0.3, max_Δt = Δt) + # initial noise to induce turbulance faster + set!(model, u = U) -simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(100)) + @info "Set initial conditions" + simulation = Simulation(model; Δt = Δt, stop_time = stop_time) -progress(sim) = @info "$(time(sim)) with Δt = $(prettytime(sim.Δt)) in $(prettytime(sim.run_wall_time))" + # Callbacks + wizard = TimeStepWizard(cfl = 0.1) + simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(100)) -simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000)) + progress(sim) = @info "$(time(sim)) with Δt = $(prettytime(sim.Δt)) in $(prettytime(sim.run_wall_time))" + simulation.callbacks[:progress] = Callback(progress, IterationInterval(1000)) -simulation.output_writers[:velocity] = JLD2OutputWriter(model, model.velocities, - overwrite_existing = true, - filename = "oscillating_cylinder_$(extra_downstream)_Re_$Re.jld2", - schedule = TimeInterval(1), - with_halos = true) + u, v, w = model.velocities + filename = "cylinder_$(simname)_Nx_$Nx.jld2" -run!(simulation) - -# load the results - -u_ts = FieldTimeSeries("oscillating_cylinder_$(extra_downstream)_Re_$Re.jld2", "u") -v_ts = FieldTimeSeries("oscillating_cylinder_$(extra_downstream)_Re_$Re.jld2", "v") - -u′, v′, w′ = Oceananigans.Fields.VelocityFields(u_ts.grid) - -ζ = Field((@at (Center, Center, Center) ∂x(v′)) - (@at (Center, Center, Center) ∂y(u′))) - -# there is probably a more memory efficient way todo this - -ζ_ts = zeros(size(grid, 1), size(grid, 2), length(u_ts.times)) # u_ts.grid so its always on cpu - -for n in 1:length(u_ts.times) - set!(u′, u_ts[n]) - set!(v′, v_ts[n]) - compute!(ζ) - ζ_ts[:, :, n] = interior(ζ, :, :, 1) + if grid isa Oceananigans.Grids.ZFlatGrid + outputs = (; model.velocities..., ζ = (@at (Center, Center, Center) ∂x(v) - ∂y(u))) + elseif grid isa Oceananigans.Grids.YFlatGrid + outputs = (; model.velocities..., ζ = (@at (Center, Center, Center) ∂x(w) - ∂z(u))) + end + simulation.output_writers[:velocity] = JLD2OutputWriter(model, outputs, + overwrite_existing = true, + filename = filename, + schedule = TimeInterval(0.5), + with_halos = true) + run!(simulation) + + if plot + # load the results + ζ_ts = FieldTimeSeries(filename, "ζ") + u_ts = FieldTimeSeries(filename, "u") + @info "Loaded results" + + xζ, yζ, zζ = nodes(ζ_ts, with_halos=true) + xu, yu, zu = nodes(u_ts, with_halos=true) + + # plot the results + fig = Figure(size = (600, 600)) + n = Observable(1) + + if grid isa Oceananigans.Grids.ZFlatGrid + ζ_plt = @lift ζ_ts[:, :, 1, $n].parent + ax = Axis(fig[1, 1], aspect = DataAspect(), xlabel = "x", ylabel = "y", width = 400, height = 400, title = "ζ") + heatmap!(ax, collect(xζ), collect(yζ), ζ_plt, colorrange = (-2, 2), colormap = :curl) + + u_plt = @lift u_ts[:, :, 1, $n].parent + ax = Axis(fig[1, 2], aspect = DataAspect(), xlabel = "x", ylabel = "y", width = 400, height = 400, title = "u") + heatmap!(ax, collect(xu), collect(yu), u_plt, colorrange = (-2, 2), colormap = :curl) + + elseif grid isa Oceananigans.Grids.YFlatGrid + ζ_plt = @lift ζ_ts[:, 1, :, $n].parent + ax = Axis(fig[1, 1], aspect = DataAspect(), xlabel = "x", ylabel = "z", width = 400, height = 400, title = "ζ") + heatmap!(ax, collect(xζ), collect(zζ), ζ_plt, colorrange = (-2, 2), colormap = :curl) + + u_plt = @lift u_ts[:, 1, :, $n].parent + ax = Axis(fig[1, 2], aspect = DataAspect(), xlabel = "x", ylabel = "z", width = 400, height = 400, title = "u") + heatmap!(ax, collect(xu), collect(zu), u_plt, colorrange = (-2, 2), colormap = :curl) + + end + resize_to_layout!(fig) + record(fig, "ζ_$filename.mp4", 1:length(ζ_ts.times), framerate = 16) do i; + n[] = i + i % 10 == 0 && @info "$(n.val) of $(length(ζ_ts.times))" + end + end end -@info "Loaded results" - -# plot the results - -fig = Figure(size = (600, 600)) - -xc, yc, zc = nodes(ζ) +matching_scheme_name(obc) = string(nameof(typeof(obc.classification.matching_scheme))) +for grid in (xygrid, xzgrid) -ax = Axis(fig[1, 1], aspect = DataAspect(), limits = (minimum(xc), maximum(xc), minimum(yc), maximum(yc))) + u_fe = FlatExtrapolationOpenBoundaryCondition(u∞, parameters = (; U, T), relaxation_timescale = 1) + v_fe = FlatExtrapolationOpenBoundaryCondition(v∞, parameters = (; U, T), relaxation_timescale = 1) + w_fe = FlatExtrapolationOpenBoundaryCondition(v∞, parameters = (; U, T), relaxation_timescale = 1) -n = Observable(1) + u_boundaries_fe = FieldBoundaryConditions(west = u_fe, east = u_fe) + v_boundaries_fe = FieldBoundaryConditions(south = v_fe, north = v_fe) + w_boundaries_fe = FieldBoundaryConditions(bottom = w_fe, top = w_fe) -ζ_plt = @lift ζ_ts[:, :, $n] - -contour!(ax, xc, yc, ζ_plt, levels = [-1, 1], colorrange = (-1, 1), colormap = :roma) + if grid isa Oceananigans.Grids.ZFlatGrid + boundary_conditions = (u = u_boundaries_fe, v = v_boundaries_fe) + simname = "xy_" * matching_scheme_name(u_boundaries_fe.east) + elseif grid isa Oceananigans.Grids.YFlatGrid + boundary_conditions = (u = u_boundaries_fe, w = w_boundaries_fe) + simname = "xz_" * matching_scheme_name(u_boundaries_fe.east) + end + run_cylinder(grid, boundary_conditions, simname = simname, stop_time = T) +end -record(fig, "oscillating_ζ_Re_$(Re)_no_exterior.mp4", 1:length(u_ts.times), framerate = 5) do i; - n[] = i - i % 10 == 0 && @info "$(n.val) of $(length(u_ts.times))" -end \ No newline at end of file