Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Default free surface for distributed grids #4061

Open
wants to merge 5 commits into
base: main
Choose a base branch
from
Open
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -49,18 +49,28 @@ mutable struct HydrostaticFreeSurfaceModel{TS, E, A<:AbstractArchitecture, S,
vertical_coordinate :: Z # Rulesets that define the time-evolution of the grid
end

default_free_surface(grid::XYRegularRG; gravitational_acceleration=g_Earth) =
default_free_surface(grid::XYRegularRG; gravitational_acceleration=g_Earth, args...) =
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you think about using the same default for all grids? I know it is maybe less efficient but there is an advantage to having the setup remain the same when you change grids / configurations. It will also probably be less confusing to debug user configurations.

When we eventually finally have tutorials on model-building, we can explain the default and how it might be optimized depending on the setting.

Basically my idea is to have a default that improves our ability to understand the model and simulations (ie maximizes convenience) --- rather than the alternative goal of choosing a default that achieves optimal performance.

Copy link
Member

@glwagner glwagner Jan 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So the first question to ask is: what is the right goal for a default?

Copy link
Collaborator Author

@simone-silvestri simone-silvestri Jan 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the right goal to strive for in a default is stability. A user wouldn't mind if the simulation is 2 or 3 times slower than it could theoretically be, it would be more annoying if he gets a NaN and doesn't know why.

So for the goal of convenience, I am ok uniformizing the free surface for all grids. We can keep the SplitExplicitFreeSurface for all grids since it is the only supported one with all grids. I don't know if we should distinguish between a variable and a fixed substep size where supported or, in the spirit of keeping all defaults the same using only a fixed substep size which is the only supported configuration for every grid. The latter option involves guessing which baroclinic time step the user will use, which might create some stability problems if the user chooses a very large time step and infringes the principle of having a stable simulation with the default configuration.

I would probably keep a variable substep size for serial grids (which should be always stable), while we could make the assumption that users that use distributed grids might have larger problems at hand and tolerate some more tampering with the settings, so as to keep the fixed substep default.

what do you think?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay, that makes sense. Basically ideally we could use the same configuration for all scenarios including Distributed, but the complexity of our performance optimization for distributed cases introduces an important trade-off. So we will compromise and use split-explicit always, but different approaches for serial vs distributed grids: variable substeps on serial grids to guarantee stability, fixed substeps on distributed grids (the only thing we support currently).

ImplicitFreeSurface(; gravitational_acceleration)

default_free_surface(grid; gravitational_acceleration=g_Earth) =
default_free_surface(grid; gravitational_acceleration=g_Earth, args...) =
SplitExplicitFreeSurface(grid; cfl = 0.7, gravitational_acceleration)

# Kind of random, probably we should infer the time-step?
default_free_surface(grid::DistributedGrid, gravitational_acceleration=g_Earth) =
SplitExplicitFreeSurface(grid; substeps = 60, gravitational_acceleration)
# A heuristic computation of a possible maximum Δt for a HydrostaticFreeSurfaceModel,
# given the grid spacings, a hypothetical `maximum_speed` (assumed to be around 3 m/s)
# and a conservative CFL of 0.3. Note that this computation only considers the ``advective''
# CFL, as we assume that at a high enough resolution advective processes are the bottleneck.
function compute_maximum_Δt(grid; maximum_speed = 3)
Δx = minimum_xspacing(grid)
Δy = minimum_yspacing(grid)
Δs = sqrt(2 / (Δx^2 + Δy^2))^(-1)
return Δs * 0.3 / maximum_speed
end

default_free_surface(grid::DistributedGrid; gravitational_acceleration=g_Earth, cfl=0.7, Δt=compute_maximum_Δt(grid)) =
SplitExplicitFreeSurface(grid; cfl=0.7, fixed_Δt=Δt, gravitational_acceleration)

default_free_surface(grid::XYRegularDistributedGrid, gravitational_acceleration=g_Earth) =
SplitExplicitFreeSurface(grid; substeps = 60, gravitational_acceleration)
default_free_surface(grid::XYRegularDistributedGrid; gravitational_acceleration=g_Earth, cfl=0.7, Δt=compute_maximum_Δt(grid)) =
SplitExplicitFreeSurface(grid; cfl=0.7, fixed_Δt=Δt, gravitational_acceleration)

"""
HydrostaticFreeSurfaceModel(; grid,
Expand All @@ -69,7 +79,10 @@ default_free_surface(grid::XYRegularDistributedGrid, gravitational_acceleration=
tracer_advection = Centered(),
buoyancy = SeawaterBuoyancy(eltype(grid)),
coriolis = nothing,
free_surface = default_free_surface(grid, gravitational_acceleration=g_Earth),
free_surface = default_free_surface(grid,
gravitational_acceleration=g_Earth,
cfl=0.7,
Δt=compute_maximum_Δt(grid)),
forcing::NamedTuple = NamedTuple(),
closure = nothing,
timestepper = :QuasiAdamsBashforth2,
Expand Down Expand Up @@ -99,7 +112,8 @@ Keyword arguments
geometry of the `grid`. If the `grid` is a `RectilinearGrid` that is
regularly spaced in the horizontal the default is an `ImplicitFreeSurface`
solver with `solver_method = :FFTBasedPoissonSolver`. In all other cases,
the default is a `SplitExplicitFreeSurface`.
the default is a `SplitExplicitFreeSurface` with a number of substeps corresponding to
a barotropic CFL of 0.7.
- `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of
preallocated `CenterField`s.
- `forcing`: `NamedTuple` of user-defined forcing functions that contribute to solution tendencies.
Expand All @@ -121,7 +135,10 @@ function HydrostaticFreeSurfaceModel(; grid,
tracer_advection = Centered(),
buoyancy = nothing,
coriolis = nothing,
free_surface = default_free_surface(grid, gravitational_acceleration=g_Earth),
free_surface = default_free_surface(grid,
gravitational_acceleration=g_Earth,
cfl=0.7,
Δt=compute_maximum_Δt(grid)),
tracers = nothing,
forcing::NamedTuple = NamedTuple(),
closure = nothing,
Expand Down
Loading