Skip to content

Commit

Permalink
Fix Corsican bug in validation test (#4070)
Browse files Browse the repository at this point in the history
* Fix a bug in validation test

* AdvectiveSkew -> AdvectiveFormulation

* Fix wall clock definition

* fix typo

* Replace WENO5 with WINO(order=5)

* an other little fix

---------

Co-authored-by: Jean-Michel Campin <[email protected]>
  • Loading branch information
glwagner and jm-c authored Jan 29, 2025
1 parent 06a7c1b commit a8cf967
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 23 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ stop_time = 30days

grid = RectilinearGrid(architecture;
topology = (Bounded, Bounded, Bounded),
size = (Ny, Ny, Nz),
size = (Ny, Ny, Nz),
x = (-Ly/2, Ly/2),
y = (-Ly/2, Ly/2),
z = (-Lz, 0),
Expand All @@ -49,8 +49,8 @@ model = HydrostaticFreeSurfaceModel(grid = grid,
buoyancy = BuoyancyTracer(),
closure = closures,
tracers = (:b, :c),
momentum_advection = WENO5(),
tracer_advection = WENO5(),
momentum_advection = WENO(order=5),
tracer_advection = WENO(order=5),
free_surface = ImplicitFreeSurface())

@info "Built $model."
Expand Down Expand Up @@ -140,7 +140,7 @@ bt = FieldTimeSeries(filepath, "b")
ct = FieldTimeSeries(filepath, "c")

# Build coordinates, rescaling the vertical coordinate
x, y, z = nodes((Center, Center, Center), grid)
x, y, z = nodes(bt)

zscale = 1
z = z .* zscale
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,11 @@ using Random
using Oceananigans
using Oceananigans.Units
using GLMakie
using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity
using Oceananigans.TurbulenceClosures: FluxTapering, DiffusiveFormulation, AdvectiveSkewClosure
using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, SkewAdvectionISSD
using Oceananigans.TurbulenceClosures: FluxTapering, DiffusiveFormulation, AdvectiveFormulation

filename = "coarse_baroclinic_adjustment"
wall_clock = Ref(time_ns())

function progress(sim)
@printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
Expand All @@ -21,7 +22,7 @@ function progress(sim)
prettytime(sim.Δt))

wall_clock[] = time_ns()

return nothing
end

Expand All @@ -37,7 +38,7 @@ save_fields_interval = 2day
stop_time = 300days

grid = RectilinearGrid(topology = (Flat, Bounded, Bounded),
size = (Ny, Nz),
size = (Ny, Nz),
y = (-Ly/2, Ly/2),
z = (-Lz, 0),
halo = (5, 5))
Expand All @@ -50,14 +51,13 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(y -> y > 0 ? -Lz : -Lz/2))
slope_limiter = FluxTapering(1000) # Allow very steep slopes

adv_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, slope_limiter)
dif_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, slope_limiter, skew_fluxes_formulation = DiffusiveFormulation())
dif_closure = IsopycnalSkewSymmetricDiffusivity(; κ_skew, slope_limiter, skew_flux_formulation = DiffusiveFormulation())

function run_simulation(closure, grid)
model = HydrostaticFreeSurfaceModel(; grid,
model = HydrostaticFreeSurfaceModel(; grid, closure,
coriolis = FPlane(latitude = -45),
buoyancy = BuoyancyTracer(),
tracer_advection = WENO(order=7),
closure = adv_closure,
tracers = (:b, :c))

@info "Built $model."
Expand Down Expand Up @@ -85,14 +85,10 @@ function run_simulation(closure, grid)
#####

simulation = Simulation(model; Δt, stop_time)
add_callback!(simulation, progress, IterationInterval(144))
suffix = closure isa SkewAdvectionISSD ? "advective" : "diffusive"

wall_clock = Ref(time_ns())

add_callback!(simulation, progress, IterationInterval(10))

suffix = closure isa AdvectiveSkewClosure ? "advective" : "diffusive"

simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
simulation.output_writers[:fields] = JLD2OutputWriter(model, merge(model.velocities, model.tracers),
schedule = TimeInterval(save_fields_interval),
filename = filename * "_fields_" * suffix,
overwrite_existing = true)
Expand Down Expand Up @@ -145,10 +141,10 @@ cd = @lift ctd[$n]

fig = Figure(size=(1800, 700))

axua = Axis(fig[2, 1], xlabel="y (km)", ylabel="z (km)", title="Zonal velocity")
axca = Axis(fig[3, 1], xlabel="y (km)", ylabel="z (km)", title="Tracer concentration")
axud = Axis(fig[2, 2], xlabel="y (km)", ylabel="z (km)", title="Zonal velocity")
axcd = Axis(fig[3, 2], xlabel="y (km)", ylabel="z (km)", title="Tracer concentration")
axua = Axis(fig[2, 1], xlabel="y (km)", ylabel="z (km)", title="GM_Adv: Zonal velocity")
axca = Axis(fig[3, 1], xlabel="y (km)", ylabel="z (km)", title="GM_Adv: Tracer concentration")
axud = Axis(fig[2, 2], xlabel="y (km)", ylabel="z (km)", title="GM_Dif: Zonal velocity")
axcd = Axis(fig[3, 2], xlabel="y (km)", ylabel="z (km)", title="GM_Dif: Tracer concentration")

levels = [-0.0015 + 0.0005 * i for i in 0:19]

Expand All @@ -174,4 +170,4 @@ display(fig)
record(fig, filename * ".mp4", 1:Nt, framerate=8) do i
@info "Plotting frame $i of $Nt"
n[] = i
end
end

0 comments on commit a8cf967

Please sign in to comment.