From 3114f74310c6de1985b7140a675ee0c79f54360e Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 14:15:23 -0800 Subject: [PATCH 01/11] this works best? --- .../one_degree_simulation/one_degree_simulation.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 08577fed..17b2dc32 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -46,9 +46,8 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) catke = ClimaOcean.OceanSimulations.default_ocean_closure() -viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000) -closure = (gm, catke, viscous_closure) +closure = (gm, catke) ##### ##### Restoring @@ -72,7 +71,9 @@ forcing = (T=FT, S=FS) ##### Ocean simulation ##### -momentum_advection = VectorInvariant() +momentum_advection = WENOVectorInvariant(vorticity_order=5, + upwinding=OnlySelfUpwinding(cross_scheme=Centered())) + tracer_advection = Centered(order=2) # Should we add a side drag since this is at a coarser resolution? @@ -95,7 +96,7 @@ atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20)) ##### coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) -simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days) +simulation = Simulation(coupled_model; Δt=2minutes, stop_time=30days) ##### ##### Run it! @@ -127,3 +128,7 @@ end add_callback!(simulation, progress, IterationInterval(10)) run!(simulation) + +simulation.Δt = 25minutes + +run!(simulation) From 9101d5791c94764cb04d477da3a20fa3ba48c6fe Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 14:18:02 -0800 Subject: [PATCH 02/11] One degree simulation take 4 --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 70c54fca..660f3321 100644 --- a/Project.toml +++ b/Project.toml @@ -42,9 +42,9 @@ JLD2 = "0.4, 0.5" KernelAbstractions = "0.9" MPI = "0.20" NCDatasets = "0.12, 0.13, 0.14" -Oceananigans = "0.94.2 - 0.99" +Oceananigans = "0.94.3 - 0.99" # This needs to be 0.94.4 when PolarBoundaryCondition is implemented OffsetArrays = "1.14" -OrthogonalSphericalShellGrids = "0.1.8" +OrthogonalSphericalShellGrids = "0.1.9" Scratch = "1" SeawaterPolynomials = "0.3.4" StaticArrays = "1" From 5833d1945fdc0a12003dd8c5cd607b6c362ef1e2 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Fri, 22 Nov 2024 15:24:30 -0800 Subject: [PATCH 03/11] add a biharmonic viscosity --- experiments/one_degree_simulation/one_degree_simulation.jl | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/experiments/one_degree_simulation/one_degree_simulation.jl b/experiments/one_degree_simulation/one_degree_simulation.jl index 17b2dc32..c551a051 100644 --- a/experiments/one_degree_simulation/one_degree_simulation.jl +++ b/experiments/one_degree_simulation/one_degree_simulation.jl @@ -46,8 +46,9 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000) catke = ClimaOcean.OceanSimulations.default_ocean_closure() +viscous_closure = HorizontalScalarBiharmonicDiffusivity(ν = 1e11) -closure = (gm, catke) +closure = (gm, catke, viscous_closure) ##### ##### Restoring @@ -71,9 +72,7 @@ forcing = (T=FT, S=FS) ##### Ocean simulation ##### -momentum_advection = WENOVectorInvariant(vorticity_order=5, - upwinding=OnlySelfUpwinding(cross_scheme=Centered())) - +momentum_advection = VectorInvariant() tracer_advection = Centered(order=2) # Should we add a side drag since this is at a coarser resolution? From 0aa74dcd863f85ff2a1b2d76158a29a7a9e0558b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Sat, 23 Nov 2024 06:27:25 -0500 Subject: [PATCH 04/11] infer maximum delta t --- src/OceanSimulations/OceanSimulations.jl | 29 +++++++++++++++++++++--- 1 file changed, 26 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 8c95965b..f87ffe4b 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -42,10 +42,33 @@ default_or_override(override, alternative_default=nothing) = override # Some defaults default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) -# 70 substeps is a safe rule of thumb for an ocean at 1/4 - 1/10th of a degree -# TODO: pass the cfl and a given Δt to calculate the number of substeps? +function infer_maximum_Δt(grid) + Δx = mean(xspacings(grid)) + Δy = mean(yspacings(grid)) + Δθ = rad2deg(mean([Δx, Δy])) / grid.radius + + # The maximum Δt is roughly 40minutes / Δθ, giving: + # - 40 minutes for a 1 degree ocean + # - 20 minutes for a 1/4 degree ocean + # - 10 minutes for a 1/8 degree ocean + # - 5 minutes for a 1/16 degree ocean + # - 2.5 minutes for a 1/32 degree ocean + + return 40minutes / Δθ +end + +# Some defaults +default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) + const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} -default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70) + +function default_free_surface(grid::TripolarOfSomeKind; + fixed_Δt = infer_maximum_Δt(grid), + cfl = 0.8) + free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) + @info "Using a $(free_surface)" + return free_surface +end function default_ocean_closure() mixing_length = CATKEMixingLength(Cᵇ=0.01) From 7c87f2b22ba43e473e13f1d228621c7eef862245 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 10:16:17 +0100 Subject: [PATCH 05/11] changes --- examples/generate_bathymetry.jl | 2 +- src/OceanSimulations/OceanSimulations.jl | 7 ++----- 2 files changed, 3 insertions(+), 6 deletions(-) diff --git a/examples/generate_bathymetry.jl b/examples/generate_bathymetry.jl index e063f379..5da3e435 100644 --- a/examples/generate_bathymetry.jl +++ b/examples/generate_bathymetry.jl @@ -38,7 +38,7 @@ grid = LatitudeLongitudeGrid(size = (Nλ, Nφ, 1), # one interpolation passes and no restrictions on connected regions. # - `h_smooth` shows the output of the function with 40 interpolation passes, which results # in a smoother bathymetry. -# - `h_no_connected_regions` shows the output of the function with `connected_regions_allowed = 0`, which +# - `h_no_connected_regions` shows the output of the function with `major_basins = 1`, which # means that the function does not allow connected regions in the bathymetry (e.g., lakes) # and fills them with land. diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index f87ffe4b..7322b375 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -57,14 +57,11 @@ function infer_maximum_Δt(grid) return 40minutes / Δθ end -# Some defaults -default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) - const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} function default_free_surface(grid::TripolarOfSomeKind; fixed_Δt = infer_maximum_Δt(grid), - cfl = 0.8) + cfl = 0.7) free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) @info "Using a $(free_surface)" return free_surface @@ -98,7 +95,7 @@ default_tracer_advection() = FluxFormAdvection(WENO(order=7), # TODO: Specify the grid to a grid on the sphere; otherwise we can provide a different # function that requires latitude and longitude etc for computing coriolis=FPlane... function ocean_simulation(grid; - Δt = 5minutes, + Δt = infer_maximum_Δt(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), From bddeae39e7fbc89305dc21c0eb5eaa308f2ce4ef Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 15:35:40 +0100 Subject: [PATCH 06/11] add statistics --- src/OceanSimulations/OceanSimulations.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 7322b375..b9098120 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -16,6 +16,7 @@ using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEEquation using SeawaterPolynomials.TEOS10: TEOS10EquationOfState +using Statistics: mean using Oceananigans.BuoyancyModels: g_Earth using Oceananigans.Coriolis: Ω_Earth From f129a1d092186696738edffefbf93c61237f53d4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 14:38:45 +0100 Subject: [PATCH 07/11] infer -> compute --- src/OceanSimulations/OceanSimulations.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 7322b375..55c67274 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -42,7 +42,7 @@ default_or_override(override, alternative_default=nothing) = override # Some defaults default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7) -function infer_maximum_Δt(grid) +function compute_maximum_Δt(grid) Δx = mean(xspacings(grid)) Δy = mean(yspacings(grid)) Δθ = rad2deg(mean([Δx, Δy])) / grid.radius @@ -60,7 +60,7 @@ end const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}} function default_free_surface(grid::TripolarOfSomeKind; - fixed_Δt = infer_maximum_Δt(grid), + fixed_Δt = compute_maximum_Δt(grid), cfl = 0.7) free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) @info "Using a $(free_surface)" @@ -95,7 +95,7 @@ default_tracer_advection() = FluxFormAdvection(WENO(order=7), # TODO: Specify the grid to a grid on the sphere; otherwise we can provide a different # function that requires latitude and longitude etc for computing coriolis=FPlane... function ocean_simulation(grid; - Δt = infer_maximum_Δt(grid), + Δt = compute_maximum_Δt(grid), closure = default_ocean_closure(), tracers = (:T, :S), free_surface = default_free_surface(grid), From c0e20384f59131167b83b56ade1ffa06d641b007 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Mon, 2 Dec 2024 14:39:15 +0100 Subject: [PATCH 08/11] compute -> infer --- src/OceanSimulations/OceanSimulations.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 55c67274..42923dec 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -16,6 +16,7 @@ using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEEquation using SeawaterPolynomials.TEOS10: TEOS10EquationOfState +using Statistics: mean using Oceananigans.BuoyancyModels: g_Earth using Oceananigans.Coriolis: Ω_Earth From 9d08000f9a52b2511bfcaab90e6de92b79d641eb Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 4 Dec 2024 18:49:07 +0100 Subject: [PATCH 09/11] default for distribtued grid --- src/OceanSimulations/OceanSimulations.jl | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 42923dec..71302d45 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -6,7 +6,7 @@ using Oceananigans using Oceananigans.Units using Oceananigans.Utils: with_tracers using Oceananigans.Advection: FluxFormAdvection -using Oceananigans.Coriolis: ActiveCellEnstrophyConserving +using Oceananigans.DistributedComputations: DistributedGrid, allreduce using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node using OrthogonalSphericalShellGrids @@ -68,6 +68,19 @@ function default_free_surface(grid::TripolarOfSomeKind; return free_surface end +function default_free_surface(grid::DistributedGrid; + fixed_Δt = compute_maximum_Δt(grid), + cfl = 0.7) + + free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt) + substeps = length(free_surface.substepping.averaging_weights) + substeps = all_reduce(max, substeps, architecture(grid)) + free_surface = SplitExplicitFreeSurface(grid; substeps) + @info "Using a $(free_surface)" + return free_surface +end + + function default_ocean_closure() mixing_length = CATKEMixingLength(Cᵇ=0.01) turbulent_kinetic_energy_equation = CATKEEquation(Cᵂϵ=1.0) From 9ead323f0e7998682650704f4c98a52b9afdfdb5 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 4 Dec 2024 19:13:59 +0100 Subject: [PATCH 10/11] all_reduce instead of allreduce --- src/OceanSimulations/OceanSimulations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/OceanSimulations/OceanSimulations.jl b/src/OceanSimulations/OceanSimulations.jl index 71302d45..7045fde1 100644 --- a/src/OceanSimulations/OceanSimulations.jl +++ b/src/OceanSimulations/OceanSimulations.jl @@ -6,7 +6,7 @@ using Oceananigans using Oceananigans.Units using Oceananigans.Utils: with_tracers using Oceananigans.Advection: FluxFormAdvection -using Oceananigans.DistributedComputations: DistributedGrid, allreduce +using Oceananigans.DistributedComputations: DistributedGrid, all_reduce using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node using OrthogonalSphericalShellGrids From 94259089f1a5b418df32399a86240ef958729c02 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 10 Dec 2024 19:52:42 +0100 Subject: [PATCH 11/11] update other packages --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index f2eec8bc..c03833c4 100644 --- a/Project.toml +++ b/Project.toml @@ -42,9 +42,9 @@ JLD2 = "0.4, 0.5" KernelAbstractions = "0.9" MPI = "0.20" NCDatasets = "0.12, 0.13, 0.14" -Oceananigans = "0.94.3 - 0.99" # This needs to be 0.94.4 when PolarBoundaryCondition is implemented +Oceananigans = "0.95 - 0.99" OffsetArrays = "1.14" -OrthogonalSphericalShellGrids = "0.1.9" +OrthogonalSphericalShellGrids = "0.1.10" Scratch = "1" SeawaterPolynomials = "0.3.4" StaticArrays = "1"