diff --git a/examples/internal_tide.jl b/examples/internal_tide.jl index 8209e80d70..44bacda8d6 100644 --- a/examples/internal_tide.jl +++ b/examples/internal_tide.jl @@ -114,12 +114,9 @@ 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. +# We built a `HydrostaticFreeSurfaceModel`: -free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7) - -model = HydrostaticFreeSurfaceModel(; grid, coriolis, free_surface, +model = HydrostaticFreeSurfaceModel(; grid, coriolis, buoyancy = BuoyancyTracer(), tracers = :b, momentum_advection = WENO(), diff --git a/src/Grids/rectilinear_grid.jl b/src/Grids/rectilinear_grid.jl index 6d22ad81e7..90deab5bbf 100644 --- a/src/Grids/rectilinear_grid.jl +++ b/src/Grids/rectilinear_grid.jl @@ -492,4 +492,3 @@ const C = Center @inline zspacings(grid::RG, ℓx, ℓy, ℓz; kwargs...) = zspacings(grid, ℓz; kwargs...) @inline isrectilinear(::RG) = true - diff --git a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl index 5c1531d8a2..74f65888e7 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl @@ -28,8 +28,8 @@ fill_horizontal_velocity_halos!(args...) = nothing ##### HydrostaticFreeSurfaceModel definition ##### -FreeSurfaceDisplacementField(velocities, free_surface, grid) = ZFaceField(grid, indices = (:, :, size(grid, 3)+1)) -FreeSurfaceDisplacementField(velocities, ::Nothing, grid) = nothing +free_surface_displacement_field(velocities, free_surface, grid) = ZFaceField(grid, indices = (:, :, size(grid, 3)+1)) +free_surface_displacement_field(velocities, ::Nothing, grid) = nothing include("compute_w_from_continuity.jl") include("rigid_lid.jl") diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_variables.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_variables.jl index f3bd5b2b27..893b70d4a5 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_variables.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_vertically_integrated_variables.jl @@ -5,7 +5,7 @@ using Oceananigans.AbstractOperations: Ax, Ay, GridMetricOperation function compute_vertically_integrated_lateral_areas!(∫ᶻ_A) # we have to account for halos when calculating Integrated areas, in case - # a periodic domain, where it is not guaranteed that ηₙ == ηₙ₊₁ + # a periodic domain, where it is not guaranteed that ηₙ == ηₙ₊₁ # 2 halos (instead of only 1) are necessary to accomodate the preconditioner field_grid = ∫ᶻ_A.xᶠᶜᶜ.grid diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_w_from_continuity.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_w_from_continuity.jl index b43920d41d..a3c590b840 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_w_from_continuity.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_w_from_continuity.jl @@ -12,7 +12,8 @@ Compute the vertical velocity ``w`` by integrating the continuity equation from w^{n+1} = -∫ [∂/∂x (u^{n+1}) + ∂/∂y (v^{n+1})] dz ``` """ -compute_w_from_continuity!(model; kwargs...) = compute_w_from_continuity!(model.velocities, model.architecture, model.grid; kwargs...) +compute_w_from_continuity!(model; kwargs...) = + compute_w_from_continuity!(model.velocities, model.architecture, model.grid; kwargs...) compute_w_from_continuity!(velocities, arch, grid; parameters = w_kernel_parameters(grid)) = launch!(arch, grid, parameters, _compute_w_from_continuity!, velocities, grid) @@ -37,11 +38,11 @@ end Hx, Hy, _ = halo_size(grid) Tx, Ty, _ = topology(grid) - Sx = Tx == Flat ? Nx : Nx + 2Hx - 2 - Sy = Ty == Flat ? Ny : Ny + 2Hy - 2 + Sx = Tx == Flat ? Nx : Nx + 2Hx - 2 + Sy = Ty == Flat ? Ny : Ny + 2Hy - 2 - Ox = Tx == Flat ? 0 : - Hx + 1 - Oy = Ty == Flat ? 0 : - Hy + 1 + Ox = Tx == Flat ? 0 : - Hx + 1 + Oy = Ty == Flat ? 0 : - Hy + 1 return KernelParameters((Sx, Sy), (Ox, Oy)) end diff --git a/src/Models/HydrostaticFreeSurfaceModels/distributed_split_explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/distributed_split_explicit_free_surface.jl index 4f14f69f2a..0dbcd5eb66 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/distributed_split_explicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/distributed_split_explicit_free_surface.jl @@ -3,7 +3,7 @@ using Oceananigans.DistributedComputations: DistributedGrid, DistributedField using Oceananigans.DistributedComputations: SynchronizedDistributed, synchronize_communication! using Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitState, SplitExplicitFreeSurface -import Oceananigans.Models.HydrostaticFreeSurfaceModels: FreeSurface, SplitExplicitAuxiliaryFields +import Oceananigans.Models.HydrostaticFreeSurfaceModels: materialize_free_surface, SplitExplicitAuxiliaryFields function SplitExplicitAuxiliaryFields(grid::DistributedGrid) @@ -59,7 +59,8 @@ end return (Ax, Ay) end -function FreeSurface(free_surface::SplitExplicitFreeSurface, velocities, grid::DistributedGrid) +# Internal function for HydrostaticFreeSurfaceModel +function materialize_free_surface(free_surface::SplitExplicitFreeSurface, velocities, grid::DistributedGrid) settings = free_surface.settings diff --git a/src/Models/HydrostaticFreeSurfaceModels/explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/explicit_free_surface.jl index 3d210e6778..5138a4f755 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/explicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/explicit_free_surface.jl @@ -25,15 +25,15 @@ Adapt.adapt_structure(to, free_surface::ExplicitFreeSurface) = ExplicitFreeSurface(Adapt.adapt(to, free_surface.η), free_surface.gravitational_acceleration) on_architecture(to, free_surface::ExplicitFreeSurface) = - ExplicitFreeSurface(on_architecture(to, free_surface.η), + ExplicitFreeSurface(on_architecture(to, free_surface.η), on_architecture(to, free_surface.gravitational_acceleration)) ##### ##### Interface to HydrostaticFreeSurfaceModel ##### -function FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, velocities, grid) - η = FreeSurfaceDisplacementField(velocities, free_surface, grid) +function materialize_free_surface(free_surface::ExplicitFreeSurface{Nothing}, velocities, grid) + η = free_surface_displacement_field(velocities, free_surface, grid) g = convert(eltype(grid), free_surface.gravitational_acceleration) return ExplicitFreeSurface(η, g) diff --git a/src/Models/HydrostaticFreeSurfaceModels/fft_based_implicit_free_surface_solver.jl b/src/Models/HydrostaticFreeSurfaceModels/fft_based_implicit_free_surface_solver.jl index 1666847009..357c258ed3 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/fft_based_implicit_free_surface_solver.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/fft_based_implicit_free_surface_solver.jl @@ -113,4 +113,3 @@ end δ_Q = flux_div_xyᶜᶜᶠ(i, j, k_top, grid, ∫ᶻQ.u, ∫ᶻQ.v) @inbounds rhs[i, j, 1] = (δ_Q - Az * η[i, j, k_top] / Δt) / (g * Lz * Δt * Az) end - diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 717bb1fd2d..a1069b231b 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -1,7 +1,6 @@ using Oceananigans.Fields: location using Oceananigans.TimeSteppers: ab2_step_field! using Oceananigans.TurbulenceClosures: implicit_step! - using Oceananigans.ImmersedBoundaries: active_interior_map, active_surface_map import Oceananigans.TimeSteppers: ab2_step! @@ -91,4 +90,3 @@ function ab2_step_tracers!(tracers, model, Δt, χ) return nothing end - diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl index 5601f85192..92f787bfd2 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_field_tuples.jl @@ -11,7 +11,7 @@ end function HydrostaticFreeSurfaceTendencyFields(velocities, free_surface, grid, tracer_names) u = XFaceField(grid) v = YFaceField(grid) - η = FreeSurfaceDisplacementField(velocities, free_surface, grid) + η = free_surface_displacement_field(velocities, free_surface, grid) tracers = TracerFields(tracer_names, grid) return merge((u=u, v=v, η=η), tracers) end diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl index b268031091..e090d20ea4 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl @@ -9,8 +9,7 @@ using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions using Oceananigans.Biogeochemistry: validate_biogeochemistry, AbstractBiogeochemistry, biogeochemical_auxiliary_fields using Oceananigans.Fields: Field, CenterField, tracernames, VelocityFields, TracerFields using Oceananigans.Forcings: model_forcing -using Oceananigans.Grids: halo_size -using Oceananigans.Grids: AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid, architecture +using Oceananigans.Grids: AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid, architecture, halo_size using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid using Oceananigans.Models: AbstractModel, validate_model_halo, NaNChecker, validate_tracer_advection, extract_boundary_conditions using Oceananigans.TimeSteppers: Clock, TimeStepper, update_state!, AbstractLagrangianParticles @@ -48,6 +47,12 @@ mutable struct HydrostaticFreeSurfaceModel{TS, E, A<:AbstractArchitecture, S, auxiliary_fields :: AF # User-specified auxiliary fields for forcing functions and boundary conditions end +default_free_surface(grid::XYRegularRG; gravitational_acceleration=g_Earth) = + ImplicitFreeSurface(; gravitational_acceleration) + +default_free_surface(grid; gravitational_acceleration=g_Earth) = + SplitExplicitFreeSurface(grid; cfl = 0.7, gravitational_acceleration) + """ HydrostaticFreeSurfaceModel(; grid, clock = Clock{eltype(grid)}(0, 0, 1), @@ -55,7 +60,7 @@ end tracer_advection = CenteredSecondOrder(), buoyancy = SeawaterBuoyancy(eltype(grid)), coriolis = nothing, - free_surface = ImplicitFreeSurface(gravitational_acceleration = g_Earth), + free_surface = default_free_surface(grid, gravitational_acceleration=g_Earth), forcing::NamedTuple = NamedTuple(), closure = nothing, boundary_conditions::NamedTuple = NamedTuple(), @@ -74,18 +79,22 @@ Keyword arguments ================= - `grid`: (required) The resolution and discrete geometry on which `model` is solved. The - architecture (CPU/GPU) that the model is solve is inferred from the architecture - of the grid. + architecture (CPU/GPU) that the model is solved is inferred from the architecture + of the `grid`. - `momentum_advection`: The scheme that advects velocities. See `Oceananigans.Advection`. - `tracer_advection`: The scheme that advects tracers. See `Oceananigans.Advection`. - `buoyancy`: The buoyancy model. See `Oceananigans.BuoyancyModels`. - `coriolis`: Parameters for the background rotation rate of the model. + - `free_surface`: The free surface model. The default free-surface solver depends on the + 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`. + - `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. - - `free_surface`: The free surface model. - `closure`: The turbulence closure for `model`. See `Oceananigans.TurbulenceClosures`. - `boundary_conditions`: `NamedTuple` containing field boundary conditions. - - `tracers`: A tuple of symbols defining the names of the modeled tracers, or a `NamedTuple` of - preallocated `CenterField`s. - `particles`: Lagrangian particles to be advected with the flow. Default: `nothing`. - `biogeochemistry`: Biogeochemical model for `tracers`. - `velocities`: The model velocities. Default: `nothing`. @@ -99,11 +108,11 @@ function HydrostaticFreeSurfaceModel(; grid, tracer_advection = CenteredSecondOrder(), buoyancy = SeawaterBuoyancy(eltype(grid)), coriolis = nothing, - free_surface = ImplicitFreeSurface(gravitational_acceleration = g_Earth), - forcing::NamedTuple = NamedTuple(), + free_surface = default_free_surface(grid, gravitational_acceleration=g_Earth), + tracers = (:T, :S), + forcing::NamedTuple = NamedTuple(), closure = nothing, boundary_conditions::NamedTuple = NamedTuple(), - tracers = (:T, :S), particles::ParticlesOrNothing = nothing, biogeochemistry::AbstractBGCOrNothing = nothing, velocities = nothing, @@ -165,7 +174,7 @@ function HydrostaticFreeSurfaceModel(; grid, @apply_regionally validate_velocity_boundary_conditions(grid, velocities) free_surface = validate_free_surface(arch, free_surface) - free_surface = FreeSurface(free_surface, velocities, grid) + free_surface = materialize_free_surface(free_surface, velocities, grid) # Instantiate timestepper if not already instantiated implicit_solver = implicit_diffusion_solver(time_discretization(closure), grid) diff --git a/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl index d745d4dc0a..ca21bfab5e 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl @@ -90,8 +90,8 @@ on_architecture(to, free_surface::ImplicitFreeSurface) = on_architecture(to, free_surface.solver_settings)) # Internal function for HydrostaticFreeSurfaceModel -function FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, grid) - η = FreeSurfaceDisplacementField(velocities, free_surface, grid) +function materialize_free_surface(free_surface::ImplicitFreeSurface{Nothing}, velocities, grid) + η = free_surface_displacement_field(velocities, free_surface, grid) gravitational_acceleration = convert(eltype(grid), free_surface.gravitational_acceleration) # Initialize barotropic volume fluxes @@ -111,13 +111,11 @@ function FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, gri free_surface.solver_settings) end -is_horizontally_regular(grid) = false -is_horizontally_regular(::RectilinearGrid{<:Any, <:Any, <:Any, <:Any, <:Number, <:Number}) = true +build_implicit_step_solver(::Val{:Default}, grid::XYRegularRG, settings, gravitational_acceleration) = + build_implicit_step_solver(Val(:FastFourierTransform), grid, settings, gravitational_acceleration) -function build_implicit_step_solver(::Val{:Default}, grid, settings, gravitational_acceleration) - default_method = is_horizontally_regular(grid) ? :FastFourierTransform : :HeptadiagonalIterativeSolver - return build_implicit_step_solver(Val(default_method), grid, settings, gravitational_acceleration) -end +build_implicit_step_solver(::Val{:Default}, grid, settings, gravitational_acceleration) = + build_implicit_step_solver(Val(:HeptadiagonalIterativeSolver), grid, settings, gravitational_acceleration) @inline explicit_barotropic_pressure_x_gradient(i, j, k, grid, ::ImplicitFreeSurface) = 0 @inline explicit_barotropic_pressure_y_gradient(i, j, k, grid, ::ImplicitFreeSurface) = 0 @@ -141,7 +139,7 @@ function implicit_free_surface_step!(free_surface::ImplicitFreeSurface, model, # Compute right hand side of implicit free surface equation @apply_regionally local_compute_integrated_volume_flux!(∫ᶻQ, model.velocities, arch) fill_halo_regions!(∫ᶻQ) - + compute_implicit_free_surface_right_hand_side!(rhs, solver, g, Δt, ∫ᶻQ, η) # Solve for the free surface at tⁿ⁺¹ @@ -157,7 +155,7 @@ function implicit_free_surface_step!(free_surface::ImplicitFreeSurface, model, end function local_compute_integrated_volume_flux!(∫ᶻQ, velocities, arch) - + foreach(mask_immersed_field!, velocities) # Compute barotropic volume flux. Blocking. diff --git a/src/Models/HydrostaticFreeSurfaceModels/matrix_implicit_free_surface_solver.jl b/src/Models/HydrostaticFreeSurfaceModels/matrix_implicit_free_surface_solver.jl index 65f4030132..ffb73f85f5 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/matrix_implicit_free_surface_solver.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/matrix_implicit_free_surface_solver.jl @@ -75,7 +75,7 @@ function solve!(η, implicit_free_surface_solver::MatrixImplicitFreeSurfaceSolve storage = implicit_free_surface_solver.storage solve!(storage, solver, rhs, Δt) - + set!(η, reshape(storage, solver.problem_size...)) return nothing @@ -92,7 +92,7 @@ function compute_implicit_free_surface_right_hand_side!(rhs, launch!(arch, grid, :xy, implicit_linearized_free_surface_right_hand_side!, rhs, grid, g, Δt, ∫ᶻQ, η) - + return nothing end @@ -130,8 +130,8 @@ end @kernel function _compute_coefficients!(diag, Ax, Ay, ∫Ax, ∫Ay, grid, g) i, j = @index(Global, NTuple) @inbounds begin - Ay[i, j, 1] = ∫Ay[i, j, 1] / Δyᶜᶠᶠ(i, j, grid.Nz+1, grid) - Ax[i, j, 1] = ∫Ax[i, j, 1] / Δxᶠᶜᶠ(i, j, grid.Nz+1, grid) - diag[i, j, 1] = - Azᶜᶜᶠ(i, j, grid.Nz+1, grid) / g + Ay[i, j, 1] = ∫Ay[i, j, 1] / Δyᶜᶠᶠ(i, j, grid.Nz+1, grid) + Ax[i, j, 1] = ∫Ax[i, j, 1] / Δxᶠᶜᶠ(i, j, grid.Nz+1, grid) + diag[i, j, 1] = - Azᶜᶜᶠ(i, j, grid.Nz+1, grid) / g end end diff --git a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl index 6d1dc41de0..471f5468f2 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/prescribed_hydrostatic_velocity_fields.jl @@ -95,12 +95,12 @@ compute_w_from_continuity!(::PrescribedVelocityFields, args...; kwargs...) = not validate_velocity_boundary_conditions(grid, ::PrescribedVelocityFields) = nothing extract_boundary_conditions(::PrescribedVelocityFields) = NamedTuple() -FreeSurfaceDisplacementField(::PrescribedVelocityFields, ::Nothing, grid) = nothing +free_surface_displacement_field(::PrescribedVelocityFields, ::Nothing, grid) = nothing HorizontalVelocityFields(::PrescribedVelocityFields, grid) = nothing, nothing -FreeSurface(::ExplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, grid) = nothing -FreeSurface(::ImplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, grid) = nothing -FreeSurface(::SplitExplicitFreeSurface, ::PrescribedVelocityFields, grid) = nothing +materialize_free_surface(::ExplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, grid) = nothing +materialize_free_surface(::ImplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, grid) = nothing +materialize_free_surface(::SplitExplicitFreeSurface, ::PrescribedVelocityFields, grid) = nothing hydrostatic_prognostic_fields(::PrescribedVelocityFields, ::Nothing, tracers) = tracers compute_hydrostatic_momentum_tendencies!(model, ::PrescribedVelocityFields, kernel_parameters; kwargs...) = nothing @@ -121,14 +121,14 @@ on_architecture(to, velocities::PrescribedVelocityFields) = # If the model only tracks particles... do nothing but that!!! const OnlyParticleTrackingModel = HydrostaticFreeSurfaceModel{TS, E, A, S, G, T, V, B, R, F, P, U, C} where - {TS, E, A, S, G, T, V, B, R, F, P<:AbstractLagrangianParticles, U<:PrescribedVelocityFields, C<:NamedTuple{(), Tuple{}}} + {TS, E, A, S, G, T, V, B, R, F, P<:AbstractLagrangianParticles, U<:PrescribedVelocityFields, C<:NamedTuple{(), Tuple{}}} -function time_step!(model::OnlyParticleTrackingModel, Δt; callbacks = [], kwargs...) +function time_step!(model::OnlyParticleTrackingModel, Δt; callbacks = [], kwargs...) model.timestepper.previous_Δt = Δt tick!(model.clock, Δt) step_lagrangian_particles!(model, Δt) update_state!(model, callbacks) end -update_state!(model::OnlyParticleTrackingModel, callbacks) = +update_state!(model::OnlyParticleTrackingModel, callbacks) = [callback(model) for callback in callbacks if callback.callsite isa UpdateStateCallsite] diff --git a/src/Models/HydrostaticFreeSurfaceModels/show_hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/show_hydrostatic_free_surface_model.jl index a423ea83f8..56b649b80b 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/show_hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/show_hydrostatic_free_surface_model.jl @@ -27,7 +27,7 @@ function Base.show(io::IO, model::HydrostaticFreeSurfaceModel) end if typeof(model.free_surface).name.wrapper == SplitExplicitFreeSurface - print(io, "│ └── substepping: $(model.free_surface.settings.substepping)", "\n") + print(io, "│ └── substepping: $(summary(model.free_surface.settings.substepping))", "\n") end end diff --git a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl index b04d251af0..e3e5d8b3c5 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl @@ -26,10 +26,12 @@ const SingleColumnGrid = AbstractGrid{<:AbstractFloat, <:Flat, <:Flat, <:Bounded ##### PressureField(arch, ::SingleColumnGrid) = (pHY′ = nothing,) -FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, velocities, ::SingleColumnGrid) = nothing -FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, ::SingleColumnGrid) = nothing -FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, ::SingleColumnGrid) = nothing -FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, ::SingleColumnGrid) = nothing +materialize_free_surface(free_surface::ExplicitFreeSurface{Nothing}, velocities, ::SingleColumnGrid) = nothing +materialize_free_surface(free_surface::ImplicitFreeSurface{Nothing}, velocities, ::SingleColumnGrid) = nothing +materialize_free_surface(free_surface::SplitExplicitFreeSurface, velocities, ::SingleColumnGrid) = nothing +materialize_free_surface(free_surface::ExplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, ::SingleColumnGrid) = nothing +materialize_free_surface(free_surface::ImplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, ::SingleColumnGrid) = nothing +materialize_free_surface(free_surface::SplitExplicitFreeSurface, ::PrescribedVelocityFields, ::SingleColumnGrid) = nothing function HydrostaticFreeSurfaceVelocityFields(::Nothing, grid::SingleColumnGrid, clock, bcs=NamedTuple()) u = XFaceField(grid, boundary_conditions=bcs.u) @@ -144,4 +146,3 @@ end @inbounds coriolis = coriolis_array[i, j] return y_f_cross_U(i, j, k, grid, coriolis, U) end - diff --git a/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl index 0f36e03fa3..5d9a00f2ed 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl @@ -27,44 +27,56 @@ struct SplitExplicitFreeSurface{𝒩, 𝒮, ℱ, 𝒫 ,ℰ} <: AbstractFreeSurfa auxiliary :: ℱ "Gravitational acceleration" gravitational_acceleration :: 𝒫 - "Settings for the split-explicit scheme (`NamedTuple`)" + "Settings for the split-explicit scheme" settings :: ℰ end """ - SplitExplicitFreeSurface(; gravitational_acceleration = g_Earth, kwargs...) + SplitExplicitFreeSurface(grid = nothing; + gravitational_acceleration = g_Earth, + substeps = nothing, + cfl = nothing, + fixed_Δt = nothing, + averaging_kernel = averaging_shape_function, + timestepper = ForwardBackwardScheme()) Return a `SplitExplicitFreeSurface` representing an explicit time discretization -of oceanic free surface dynamics with `gravitational_acceleration`. +of a free surface dynamics with `gravitational_acceleration`. Keyword Arguments ================= +- `gravitational_acceleration`: the gravitational acceleration (default: `g_Earth`) + - `substeps`: The number of substeps that divide the range `(t, t + 2Δt)`, where `Δt` is the baroclinic timestep. Note that some averaging functions do not require substepping until `2Δt`. - The number of substeps is reduced automatically to the last index of `averaging_weights` - for which `averaging_weights > 0`. + The number of substeps is reduced automatically to the last index of `averaging_kernel` + for which `averaging_kernel > 0`. -- `cfl`: If set then the number of `substeps` are computed based on the advective timescale imposed from the - barotropic gravity-wave speed, computed with depth `grid.Lz`. If `fixed_Δt` is provided then the number of - `substeps` will adapt to maintain an exact cfl. If not the effective cfl will be always lower than the - specified `cfl` provided that the baroclinic time step `Δt_baroclinic < fixed_Δt` +- `cfl`: If set then the number of `substeps` are computed based on the advective timescale imposed from + the barotropic gravity-wave speed that corresponds to depth `grid.Lz`. If `fixed_Δt` is provided, + then the number of `substeps` adapts to maintain an exact `cfl`. If not, the effective cfl will + always be lower than the specified `cfl` provided that the baroclinic time step satisfies + `Δt_baroclinic < fixed_Δt`. !!! info "Needed keyword arguments" - Either `substeps` _or_ `cfl` (with `grid`) need to be prescribed. - -- `grid`: Used to compute the corresponding barotropic surface wave speed. - -- `fixed_Δt`: The maximum baroclinic timestep allowed. If `fixed_Δt` is a `nothing` and a cfl is provided, then - the number of substeps will be computed on the fly from the baroclinic time step to maintain a constant cfl. - -- `gravitational_acceleration`: the gravitational acceleration (default: `g_Earth`) - -- `averaging_kernel`: function of `τ` used to average the barotropic transport `U` and free surface `η` - within the barotropic advancement. `τ` is the fractional substep going from 0 to 2 - with the baroclinic time step `t + Δt` located at `τ = 1`. This function should be - centered at `τ = 1`, that is, ``∑ (aₘ m /M) = 1``. By default the averaging kernel - described by [Shchepetkin2005](@citet) is chosen. + Either `substeps` _or_ `cfl` need to be prescribed. + + When `cfl` is prescribed then `grid` is also required as a positional argument. + +- `fixed_Δt`: The maximum baroclinic timestep allowed. If `fixed_Δt` is a `nothing` and a cfl is provided, + then the number of substeps will be computed on the fly from the baroclinic time step to + maintain a constant cfl. + +- `averaging_kernel`: A function of `τ` used to average the barotropic transport `U` and the free surface + `η` within the barotropic advancement. `τ` is the fractional substep going from 0 to 2 + with the baroclinic time step `t + Δt` located at `τ = 1`. The `averaging_kernel` + function should be centered at `τ = 1`, that is, ``∑ (aₘ m / M) = 1``, where the + the summation occurs for ``m = 1, ..., M_*``. Here, ``m = 0`` and ``m = M`` correspond + to the two consecutive baroclinic timesteps between which the barotropic timestepping + occurs and ``M_*`` corresponds to the last barotropic time step for which the + `averaging_kernel > 0`. By default, the averaging kernel described by [Shchepetkin2005](@citet) + is used. - `timestepper`: Time stepping scheme used for the barotropic advancement. Choose one of: * `ForwardBackwardScheme()` (default): `η = f(U)` then `U = f(η)`, @@ -75,31 +87,45 @@ References Shchepetkin, A. F., & McWilliams, J. C. (2005). The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modelling, 9(4), 347-404. """ -SplitExplicitFreeSurface(FT::DataType = Float64; gravitational_acceleration = g_Earth, kwargs...) = - SplitExplicitFreeSurface(nothing, nothing, nothing, convert(FT, gravitational_acceleration), - SplitExplicitSettings(FT; gravitational_acceleration, kwargs...)) - -# The new constructor is defined later on after the state, settings, auxiliary have been defined -function FreeSurface(free_surface::SplitExplicitFreeSurface, velocities, grid) - η = FreeSurfaceDisplacementField(velocities, free_surface, grid) - - return SplitExplicitFreeSurface(η, SplitExplicitState(grid, free_surface.settings.timestepper), - SplitExplicitAuxiliaryFields(grid), - free_surface.gravitational_acceleration, - free_surface.settings) +function SplitExplicitFreeSurface(grid = nothing; + gravitational_acceleration = g_Earth, + substeps = nothing, + cfl = nothing, + fixed_Δt = nothing, + averaging_kernel = averaging_shape_function, + timestepper = ForwardBackwardScheme()) + + settings = SplitExplicitSettings(grid; + gravitational_acceleration, + substeps, + cfl, + fixed_Δt, + averaging_kernel, + timestepper) + + return SplitExplicitFreeSurface(nothing, + nothing, + nothing, + gravitational_acceleration, + settings) end -function SplitExplicitFreeSurface(grid; gravitational_acceleration = g_Earth, - settings = SplitExplicitSettings(eltype(grid); gravitational_acceleration, substeps = 200)) +# Internal function for HydrostaticFreeSurfaceModel +function materialize_free_surface(free_surface::SplitExplicitFreeSurface, velocities, grid) + settings = SplitExplicitSettings(grid; free_surface.settings.settings_kwargs...) - Nz = size(grid, 3) - η = ZFaceField(grid, indices = (:, :, Nz+1)) - gravitational_acceleration = convert(eltype(grid), gravitational_acceleration) + η = free_surface_displacement_field(velocities, free_surface, grid) - return SplitExplicitFreeSurface(η, SplitExplicitState(grid, settings.timestepper), SplitExplicitAuxiliaryFields(grid), - gravitational_acceleration, settings) + gravitational_acceleration = convert(eltype(grid), free_surface.gravitational_acceleration) + + return SplitExplicitFreeSurface(η, + SplitExplicitState(grid, settings.timestepper), + SplitExplicitAuxiliaryFields(grid), + gravitational_acceleration, + settings) end + """ struct SplitExplicitState @@ -144,15 +170,15 @@ acts as a filter for `η`. Values with superscripts `m-1` and `m-2` correspond t time steps to allow using a higher-order time stepping scheme, e.g., `AdamsBashforth3Scheme`. """ function SplitExplicitState(grid::AbstractGrid, timestepper) - + Nz = size(grid, 3) - + η̅ = ZFaceField(grid, indices = (:, :, Nz+1)) ηᵐ = auxiliary_free_surface_field(grid, timestepper) ηᵐ⁻¹ = auxiliary_free_surface_field(grid, timestepper) ηᵐ⁻² = auxiliary_free_surface_field(grid, timestepper) - + U = XFaceField(grid, indices = (:, :, Nz)) V = YFaceField(grid, indices = (:, :, Nz)) @@ -160,15 +186,15 @@ function SplitExplicitState(grid::AbstractGrid, timestepper) Vᵐ⁻¹ = auxiliary_barotropic_V_field(grid, timestepper) Uᵐ⁻² = auxiliary_barotropic_U_field(grid, timestepper) Vᵐ⁻² = auxiliary_barotropic_V_field(grid, timestepper) - + U̅ = XFaceField(grid, indices = (:, :, Nz)) V̅ = YFaceField(grid, indices = (:, :, Nz)) - + return SplitExplicitState(; ηᵐ, ηᵐ⁻¹, ηᵐ⁻², U, Uᵐ⁻¹, Uᵐ⁻², V, Vᵐ⁻¹, Vᵐ⁻², η̅, U̅, V̅) end """ - SplitExplicitAuxiliaryFields + struct SplitExplicitAuxiliaryFields A type containing auxiliary fields for the split-explicit free surface. @@ -225,8 +251,9 @@ A type containing settings for the split-explicit free surface. $(FIELDS) """ struct SplitExplicitSettings{𝒩, 𝒮} - substepping :: 𝒩 # Either `FixedSubstepNumber` or `FixedTimeStepSize`" - timestepper :: 𝒮 # time-stepping scheme + substepping :: 𝒩 # Either `FixedSubstepNumber` or `FixedTimeStepSize`" + timestepper :: 𝒮 # time-stepping scheme + settings_kwargs :: NamedTuple # kwargs to reproduce current settings end struct AdamsBashforth3Scheme end @@ -242,7 +269,7 @@ auxiliary_barotropic_V_field(grid, ::AdamsBashforth3Scheme) = YFaceField(grid, i auxiliary_barotropic_V_field(grid, ::ForwardBackwardScheme) = nothing # (p = 2, q = 4, r = 0.18927) minimize dispersion error from Shchepetkin and McWilliams (2005): https://doi.org/10.1016/j.ocemod.2004.08.002 -@inline function averaging_shape_function(τ::FT; p = 2, q = 4, r = FT(0.18927)) where FT +@inline function averaging_shape_function(τ::FT; p = 2, q = 4, r = FT(0.18927)) where FT τ₀ = (p + 2) * (p + q + 2) / (p + 1) / (p + q + 1) return (τ / τ₀)^p * (1 - (τ / τ₀)^q) - r * (τ / τ₀) @@ -252,7 +279,7 @@ end @inline constant_averaging_kernel(τ::FT) where FT = convert(FT, 1) """ An internal type for the `SplitExplicitFreeSurface` that allows substepping with -a fixed `Δt_barotopic` based on a CFL condition """ +a fixed `Δt_barotropic` based on a CFL condition """ struct FixedTimeStepSize{B, F} Δt_barotropic :: B averaging_kernel :: F @@ -265,11 +292,13 @@ struct FixedSubstepNumber{B, F} averaging_weights :: F end -function FixedTimeStepSize(FT::DataType = Float64; - cfl = 0.7, - grid, - averaging_kernel = averaging_shape_function, +function FixedTimeStepSize(grid; + cfl = 0.7, + averaging_kernel = averaging_shape_function, gravitational_acceleration = g_Earth) + + FT = eltype(grid) + Δx⁻² = topology(grid)[1] == Flat ? 0 : 1 / minimum_xspacing(grid)^2 Δy⁻² = topology(grid)[2] == Flat ? 0 : 1 / minimum_yspacing(grid)^2 Δs = sqrt(1 / (Δx⁻² + Δy⁻²)) @@ -296,30 +325,43 @@ end return Δτ, tuple(averaging_weights...) end -function SplitExplicitSettings(FT::DataType=Float64; - substeps = nothing, - cfl = nothing, - grid = nothing, - fixed_Δt = nothing, +function SplitExplicitSettings(grid = nothing; gravitational_acceleration = g_Earth, + substeps = nothing, + cfl = nothing, + fixed_Δt = nothing, averaging_kernel = averaging_shape_function, timestepper = ForwardBackwardScheme()) - - if (!isnothing(substeps) && !isnothing(cfl)) || (isnothing(substeps) && isnothing(cfl)) - throw(ArgumentError("either specify a cfl or a number of substeps")) + + settings_kwargs = (; gravitational_acceleration, + substeps, + cfl, + fixed_Δt, + averaging_kernel, + timestepper) + + if !isnothing(grid) + FT = eltype(grid) + else + # this is a fallback and only used via the outer constructor, + # in case no grid is provided; when afterwards the free surfade + # is materialized via materialize_free_surface + # FT becomes eltype(grid) + FT = Float64 end - if !isnothing(grid) && eltype(grid) !== FT - throw(ArgumentError("Prescribed FT was different that the one used in `grid`.")) + if (!isnothing(substeps) && !isnothing(cfl)) || (isnothing(substeps) && isnothing(cfl)) + throw(ArgumentError("either specify a cfl or a number of substeps")) end if !isnothing(cfl) if isnothing(grid) - throw(ArgumentError("Need to specify the grid kwarg to calculate the barotropic substeps from the cfl")) + throw(ArgumentError(string("Need to provide the grid to calculate the barotropic substeps from the cfl. ", + "For example, SplitExplicitFreeSurface(grid, cfl=0.7, ...)"))) end - substepping = FixedTimeStepSize(FT; cfl, grid, gravitational_acceleration, averaging_kernel) + substepping = FixedTimeStepSize(grid; cfl, gravitational_acceleration, averaging_kernel) if isnothing(fixed_Δt) - return SplitExplicitSettings(substepping, timestepper) + return SplitExplicitSettings(substepping, timestepper, settings_kwargs) else substeps = ceil(Int, 2 * fixed_Δt / substepping.Δt_barotropic) end @@ -328,7 +370,7 @@ function SplitExplicitSettings(FT::DataType=Float64; fractional_step_size, averaging_weights = weights_from_substeps(FT, substeps, averaging_kernel) substepping = FixedSubstepNumber(fractional_step_size, averaging_weights) - return SplitExplicitSettings(substepping, timestepper) + return SplitExplicitSettings(substepping, timestepper, settings_kwargs) end # Convenience Functions for grabbing free surface @@ -342,10 +384,11 @@ free_surface(free_surface::SplitExplicitFreeSurface) = free_surface.η (sefs::SplitExplicitFreeSurface)(settings::SplitExplicitSettings) = SplitExplicitFreeSurface(sefs.η, sefs.state, sefs.auxiliary, sefs.gravitational_acceleration, settings) -Base.summary(s::FixedTimeStepSize) = string("Barotropic time step equal to $(s.Δt_barotopic)") +Base.summary(s::FixedTimeStepSize) = string("Barotropic time step equal to $(prettytime(s.Δt_barotropic))") Base.summary(s::FixedSubstepNumber) = string("Barotropic fractional step equal to $(s.fractional_step_size) times the baroclinic step") -Base.summary(sefs::SplitExplicitFreeSurface) = string("SplitExplicitFreeSurface with $(sefs.settings.substepping)") +Base.summary(sefs::SplitExplicitFreeSurface) = string("SplitExplicitFreeSurface with $(summary(sefs.settings.substepping))") + Base.show(io::IO, sefs::SplitExplicitFreeSurface) = print(io, "$(summary(sefs))\n") function reset!(sefs::SplitExplicitFreeSurface) @@ -365,19 +408,17 @@ Adapt.adapt_structure(to, free_surface::SplitExplicitFreeSurface) = SplitExplicitFreeSurface(Adapt.adapt(to, free_surface.η), nothing, nothing, free_surface.gravitational_acceleration, nothing) -for Type in (:SplitExplicitFreeSurface, - :SplitExplicitSettings, +for Type in (:SplitExplicitFreeSurface, + :SplitExplicitSettings, :SplitExplicitState, :SplitExplicitAuxiliaryFields, :FixedTimeStepSize, :FixedSubstepNumber) @eval begin - function on_architecture(to, settings::$Type) - args = Tuple(on_architecture(to, prop) for prop in propertynames(settings)) - return SplitExplicitState(args...) + function on_architecture(to, fs::$Type) + args = Tuple(on_architecture(to, prop) for prop in propertynames(fs)) + return $Type(args...) end end end - - \ No newline at end of file diff --git a/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl b/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl index da8839732d..1700126983 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface_kernels.jl @@ -365,7 +365,7 @@ function iterate_split_explicit!(free_surface, grid, Δτᴮ, weights, ::Val{Nsu auxiliary = free_surface.auxiliary settings = free_surface.settings g = free_surface.gravitational_acceleration - + # unpack state quantities, parameters and forcing terms U, V = state.U, state.V Uᵐ⁻¹, Uᵐ⁻² = state.Uᵐ⁻¹, state.Uᵐ⁻² @@ -373,11 +373,11 @@ function iterate_split_explicit!(free_surface, grid, Δτᴮ, weights, ::Val{Nsu ηᵐ, ηᵐ⁻¹, ηᵐ⁻² = state.ηᵐ, state.ηᵐ⁻¹, state.ηᵐ⁻² η̅, U̅, V̅ = state.η̅, state.U̅, state.V̅ Gᵁ, Gⱽ, Hᶠᶜ, Hᶜᶠ = auxiliary.Gᵁ, auxiliary.Gⱽ, auxiliary.Hᶠᶜ, auxiliary.Hᶜᶠ - - timestepper = settings.timestepper - + + timestepper = settings.timestepper + parameters = auxiliary.kernel_parameters - + free_surface_kernel! = configured_kernel(arch, grid, parameters, _split_explicit_free_surface!) barotropic_velocity_kernel! = configured_kernel(arch, grid, parameters, _split_explicit_barotropic_velocity!) @@ -398,7 +398,7 @@ function iterate_split_explicit!(free_surface, grid, Δτᴮ, weights, ::Val{Nsu # To alleviate this penalty we convert first and then we substep! converted_η_args = convert_args(arch, η_args) converted_U_args = convert_args(arch, U_args) - + @unroll for substep in 1:Nsubsteps Base.@_inline_meta averaging_weight = weights[substep] @@ -410,9 +410,9 @@ function iterate_split_explicit!(free_surface, grid, Δτᴮ, weights, ::Val{Nsu return nothing end -# Calculate RHS for the barotopic time step. +# Calculate RHS for the barotropic time step. @kernel function _compute_integrated_ab2_tendencies!(Gᵁ, Gⱽ, grid, Gu⁻, Gv⁻, Guⁿ, Gvⁿ, χ) - i, j = @index(Global, NTuple) + i, j = @index(Global, NTuple) k_top = grid.Nz + 1 @inbounds Gᵁ[i, j, k_top-1] = Δzᶠᶜᶜ(i, j, 1, grid) * ab2_step_Gu(i, j, 1, grid, Gu⁻, Guⁿ, χ) @@ -424,7 +424,7 @@ end end end -# Calculate RHS for the barotopic time step. +# Calculate RHS for the barotropic time step.q @kernel function _compute_integrated_ab2_tendencies!(Gᵁ, Gⱽ, grid::ActiveZColumnsIBG, Gu⁻, Gv⁻, Guⁿ, Gvⁿ, χ) idx = @index(Global, Linear) i, j = active_linear_index_to_tuple(idx, ZColumnMap(), grid) @@ -445,7 +445,7 @@ end @inline ab2_step_Gv(i, j, k, grid, G⁻, Gⁿ, χ::FT) where FT = @inbounds ifelse(peripheral_node(i, j, k, grid, c, f, c), zero(grid), (convert(FT, 1.5) + χ) * Gⁿ[i, j, k] - G⁻[i, j, k] * (convert(FT, 0.5) + χ)) -# Setting up the RHS for the barotropic step (tendencies of the barotopic velocity components) +# Setting up the RHS for the barotropic step (tendencies of the barotropic velocity components) # This function is called after `calculate_tendency` and before `ab2_step_velocities!` function setup_free_surface!(model, free_surface::SplitExplicitFreeSurface, χ) diff --git a/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl b/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl index 4cf7c19521..c1e477138e 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl @@ -31,7 +31,7 @@ update_state!(model::HydrostaticFreeSurfaceModel, callbacks=[]; compute_tendenci function update_state!(model::HydrostaticFreeSurfaceModel, grid, callbacks; compute_tendencies = true) @apply_regionally mask_immersed_model_fields!(model, grid) - + # Update possible FieldTimeSeries used in the model @apply_regionally update_model_field_time_series!(model, model.clock) @@ -42,10 +42,10 @@ function update_state!(model::HydrostaticFreeSurfaceModel, grid, callbacks; comp fill_halo_regions!(model.diffusivity_fields; only_local_halos = true) [callback(model) for callback in callbacks if callback.callsite isa UpdateStateCallsite] - + update_biogeochemical_state!(model.biogeochemistry, model) - compute_tendencies && + compute_tendencies && @apply_regionally compute_tendencies!(model, callbacks) return nothing @@ -68,8 +68,8 @@ end function compute_auxiliaries!(model::HydrostaticFreeSurfaceModel; w_parameters = tuple(w_kernel_parameters(model.grid)), p_parameters = tuple(p_kernel_parameters(model.grid)), - κ_parameters = tuple(:xyz)) - + κ_parameters = tuple(:xyz)) + grid = model.grid closure = model.closure diffusivity = model.diffusivity_fields @@ -77,9 +77,9 @@ function compute_auxiliaries!(model::HydrostaticFreeSurfaceModel; w_parameters = for (wpar, ppar, κpar) in zip(w_parameters, p_parameters, κ_parameters) compute_w_from_continuity!(model; parameters = wpar) compute_diffusivities!(diffusivity, closure, model; parameters = κpar) - update_hydrostatic_pressure!(model.pressure.pHY′, architecture(grid), - grid, model.buoyancy, model.tracers; - parameters = ppar) + update_hydrostatic_pressure!(model.pressure.pHY′, architecture(grid), + grid, model.buoyancy, model.tracers; + parameters = ppar) end return nothing end diff --git a/src/Models/Models.jl b/src/Models/Models.jl index f7c871287b..fa46608b96 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -10,19 +10,18 @@ export seawater_density using Oceananigans: AbstractModel, fields, prognostic_fields +using Oceananigans.AbstractOperations: AbstractOperation +using Oceananigans.Advection: AbstractAdvectionScheme, CenteredSecondOrder, VectorInvariant +using Oceananigans.Fields: AbstractField, Field, flattened_unique_values, boundary_conditions using Oceananigans.Grids: AbstractGrid, halo_size, inflate_halo_size +using Oceananigans.OutputReaders: update_field_time_series!, extract_field_time_series using Oceananigans.TimeSteppers: AbstractTimeStepper, Clock using Oceananigans.Utils: Time -using Oceananigans.Fields: AbstractField, Field, flattened_unique_values, boundary_conditions -using Oceananigans.AbstractOperations: AbstractOperation -using Oceananigans.Advection: AbstractAdvectionScheme, CenteredSecondOrder, VectorInvariant import Oceananigans: initialize! import Oceananigans.Architectures: architecture import Oceananigans.TimeSteppers: reset! -using Oceananigans.OutputReaders: update_field_time_series!, extract_field_time_series - # A prototype interface for AbstractModel. # # TODO: decide if we like this. diff --git a/src/MultiRegion/multi_region_grid.jl b/src/MultiRegion/multi_region_grid.jl index 707de5bd3c..225134031d 100644 --- a/src/MultiRegion/multi_region_grid.jl +++ b/src/MultiRegion/multi_region_grid.jl @@ -3,6 +3,7 @@ using Oceananigans.ImmersedBoundaries: GridFittedBottom, PartialCellBottom, Grid import Oceananigans.Grids: architecture, size, new_data, halo_size import Oceananigans.Grids: with_halo, on_architecture +import Oceananigans.Models.HydrostaticFreeSurfaceModels: default_free_surface import Oceananigans.DistributedComputations: reconstruct_global_grid import Oceananigans.Grids: minimum_spacing, destantiate @@ -48,6 +49,10 @@ minimum_spacing(dir, grid::MultiRegionGrid, ℓx, ℓy, ℓz) = @inline Base.length(mrg::MultiRegionGrid) = Base.length(mrg.region_grids) @inline Base.length(mrg::ImmersedMultiRegionGrid) = Base.length(mrg.underlying_grid.region_grids) +# the default free surface solver; see Models.HydrostaticFreeSurfaceModels +default_free_surface(grid::MultiRegionGrid; gravitational_acceleration=g_Earth) = + ImplicitFreeSurface(; gravitational_acceleration) + """ MultiRegionGrid(global_grid; partition = XPartition(2), devices = nothing, diff --git a/src/MultiRegion/multi_region_split_explicit_free_surface.jl b/src/MultiRegion/multi_region_split_explicit_free_surface.jl index d33f6bd1a6..67246b3615 100644 --- a/src/MultiRegion/multi_region_split_explicit_free_surface.jl +++ b/src/MultiRegion/multi_region_split_explicit_free_surface.jl @@ -1,11 +1,13 @@ using Oceananigans.Utils using Oceananigans.AbstractOperations: GridMetricOperation, Δz -using Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitState, SplitExplicitFreeSurface +using Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitFreeSurface, + SplitExplicitSettings, + SplitExplicitState -import Oceananigans.Models.HydrostaticFreeSurfaceModels: FreeSurface, SplitExplicitAuxiliaryFields +import Oceananigans.Models.HydrostaticFreeSurfaceModels: materialize_free_surface, SplitExplicitAuxiliaryFields function SplitExplicitAuxiliaryFields(grid::MultiRegionGrids) - + Gᵁ = Field((Face, Center, Nothing), grid) Gⱽ = Field((Center, Face, Nothing), grid) @@ -20,7 +22,7 @@ function SplitExplicitAuxiliaryFields(grid::MultiRegionGrids) # In a non-parallel grid we calculate only the interior @apply_regionally kernel_size = augmented_kernel_size(grid, grid.partition) @apply_regionally kernel_offsets = augmented_kernel_offsets(grid, grid.partition) - + @apply_regionally kernel_parameters = KernelParameters(kernel_size, kernel_offsets) return SplitExplicitAuxiliaryFields(Gᵁ, Gⱽ, Hᶠᶜ, Hᶜᶠ, kernel_parameters) @@ -40,7 +42,10 @@ end @inline augmented_kernel_offsets(grid, ::YPartition) = (0, - halo_size(grid)[2] + 1) @inline augmented_kernel_offsets(grid, ::CubedSpherePartition) = (- halo_size(grid)[2] + 1, - halo_size(grid)[2] + 1) -function FreeSurface(free_surface::SplitExplicitFreeSurface, velocities, grid::MultiRegionGrids) +# Internal function for HydrostaticFreeSurfaceModel +function materialize_free_surface(free_surface::SplitExplicitFreeSurface, velocities, grid::MultiRegionGrids) + + settings = SplitExplicitSettings(grid; free_surface.settings.settings_kwargs...) switch_device!(grid.devices[1]) old_halos = halo_size(getregion(grid, 1)) diff --git a/src/Operators/spacings_and_areas_and_volumes.jl b/src/Operators/spacings_and_areas_and_volumes.jl index 35ef0471e4..5735444c63 100644 --- a/src/Operators/spacings_and_areas_and_volumes.jl +++ b/src/Operators/spacings_and_areas_and_volumes.jl @@ -1,9 +1,9 @@ using Oceananigans.Grids: Center, Face -const RG = RectilinearGrid -const RGX = XRegularRG -const RGY = YRegularRG -const RGZ = ZRegularRG +const RG = RectilinearGrid +const RGX = XRegularRG +const RGY = YRegularRG +const RGZ = ZRegularRG const OSSG = OrthogonalSphericalShellGrid const OSSGZ = ZRegOrthogonalSphericalShellGrid @@ -53,10 +53,10 @@ The operators in this file fall into three categories: ##### ##### -@inline Δxᶠᵃᵃ(i, j, k, grid) = nothing -@inline Δxᶜᵃᵃ(i, j, k, grid) = nothing -@inline Δyᵃᶠᵃ(i, j, k, grid) = nothing -@inline Δyᵃᶜᵃ(i, j, k, grid) = nothing +@inline Δxᶠᵃᵃ(i, j, k, grid) = nothing +@inline Δxᶜᵃᵃ(i, j, k, grid) = nothing +@inline Δyᵃᶠᵃ(i, j, k, grid) = nothing +@inline Δyᵃᶜᵃ(i, j, k, grid) = nothing const ZRG = Union{LLGZ, RGZ} diff --git a/src/Solvers/fft_based_poisson_solver.jl b/src/Solvers/fft_based_poisson_solver.jl index 22c1eddce0..084562471c 100644 --- a/src/Solvers/fft_based_poisson_solver.jl +++ b/src/Solvers/fft_based_poisson_solver.jl @@ -74,7 +74,7 @@ end """ solve!(ϕ, solver::FFTBasedPoissonSolver, b, m=0) -Solves the "generalized" Poisson equation, +Solve the "generalized" Poisson equation, ```math (∇² + m) ϕ = b, diff --git a/test/test_distributed_hydrostatic_model.jl b/test/test_distributed_hydrostatic_model.jl index 7c75436d3f..52ac851ebc 100644 --- a/test/test_distributed_hydrostatic_model.jl +++ b/test/test_distributed_hydrostatic_model.jl @@ -38,7 +38,7 @@ end function solid_body_rotation_test(grid) - free_surface = SplitExplicitFreeSurface(; substeps = 5, gravitational_acceleration = 1) + free_surface = SplitExplicitFreeSurface(grid; substeps = 5, gravitational_acceleration = 1) coriolis = HydrostaticSphericalCoriolis(rotation_rate = 1) model = HydrostaticFreeSurfaceModel(; grid, diff --git a/test/test_hydrostatic_regression.jl b/test/test_hydrostatic_regression.jl index d555ec4d81..f7937d2eaa 100644 --- a/test/test_hydrostatic_regression.jl +++ b/test/test_hydrostatic_regression.jl @@ -8,10 +8,10 @@ function show_hydrostatic_test(grid, free_surface, precompute_metrics) typeof(grid) <: XRegularLLG ? gx = :regular : gx = :stretched typeof(grid) <: YRegularLLG ? gy = :regular : gy = :stretched typeof(grid) <: ZRegularLLG ? gz = :regular : gz = :stretched - + arch = grid.architecture free_surface_str = string(typeof(free_surface).name.wrapper) - + strc = "$(precompute_metrics ? ", metrics are precomputed" : "")" testset_str = "Hydrostatic free turbulence regression [$(arch), $(topology(grid, 1)) longitude, ($gx, $gy, $gz) grid, $free_surface_str]" * strc @@ -65,27 +65,29 @@ include("regression_tests/hydrostatic_free_turbulence_regression_test.jl") longitudes = [(-180, 180), (-160, 160)] latitudes = [(-60, 60)] zs = [(-90, 0)] - + explicit_free_surface = ExplicitFreeSurface(gravitational_acceleration = 1.0) implicit_free_surface = ImplicitFreeSurface(gravitational_acceleration = 1.0, solver_method = :PreconditionedConjugateGradient, reltol = 0, abstol = 1e-15) - split_explicit_free_surface = SplitExplicitFreeSurface(gravitational_acceleration = 1.0, substeps = 5) - for longitude in longitudes, latitude in latitudes, z in zs, precompute_metrics in (true, false) longitude[1] == -180 ? size = (180, 60, 3) : size = (160, 60, 3) grid = LatitudeLongitudeGrid(arch; size, longitude, latitude, z, precompute_metrics, halo=(2, 2, 2)) + split_explicit_free_surface = SplitExplicitFreeSurface(grid, + gravitational_acceleration = 1.0, + substeps = 5) + for free_surface in [explicit_free_surface, implicit_free_surface, split_explicit_free_surface] - + # GPU + ImplicitFreeSurface + precompute metrics cannot be tested on sverdrup at the moment # because "uses too much parameter space (maximum 0x1100 bytes)" error if !(precompute_metrics && free_surface isa ImplicitFreeSurface && arch isa GPU) && !(free_surface isa ImplicitFreeSurface && arch isa Distributed) # Also no implicit free surface on distributed testset_str, info_str = show_hydrostatic_test(grid, free_surface, precompute_metrics) - + @testset "$testset_str" begin @info "$info_str" run_hydrostatic_free_turbulence_regression_test(grid, free_surface) diff --git a/test/test_immersed_advection.jl b/test/test_immersed_advection.jl index 80938b7ef2..fa948add2b 100644 --- a/test/test_immersed_advection.jl +++ b/test/test_immersed_advection.jl @@ -36,9 +36,9 @@ end function run_tracer_conservation_test(grid, scheme) - model = HydrostaticFreeSurfaceModel(grid = grid, tracers = :c, + model = HydrostaticFreeSurfaceModel(grid = grid, tracers = :c, free_surface = ExplicitFreeSurface(), - tracer_advection = scheme, + tracer_advection = scheme, buoyancy = nothing, coriolis = nothing) @@ -55,7 +55,7 @@ function run_tracer_conservation_test(grid, scheme) wave_speed = sqrt(model.free_surface.gravitational_acceleration) dt = 0.1 / wave_speed - for i in 1:10 + for _ in 1:10 time_step!(model, dt) end @@ -100,15 +100,15 @@ for arch in archs grid = RectilinearGrid(arch, size=(20, 20), extent=(20, 20), halo = (6, 6), topology=(Bounded, Bounded, Flat)) ibg = ImmersedBoundaryGrid(grid, GridFittedBoundary((x, y) -> (x < 5 || y < 5))) - + c = CenterField(ibg) set!(c, 1) mask_immersed_field!(c) fill_halo_regions!(c) - + for adv in advection_schemes, buffer in [1, 2, 3, 4, 5] scheme = adv(order = advective_order(buffer, adv)) - + @info " Testing immersed tracer reconstruction [$(typeof(arch)), $(summary(scheme))]" run_tracer_interpolation_test(c, ibg, scheme) end @@ -149,10 +149,9 @@ for arch in archs for adv in advection_schemes, buffer in [1, 2, 3, 4, 5] scheme = adv(order = advective_order(buffer, adv)) - + @info " Testing immersed momentum reconstruction [$(typeof(arch)), $(summary(scheme))]" run_momentum_interpolation_test(u, v, ibg, scheme) end end end - diff --git a/test/test_implicit_free_surface_solver.jl b/test/test_implicit_free_surface_solver.jl index e5f32ff3af..2e4bc55a2f 100644 --- a/test/test_implicit_free_surface_solver.jl +++ b/test/test_implicit_free_surface_solver.jl @@ -6,7 +6,6 @@ using Oceananigans.Operators using Oceananigans.Grids: inactive_cell using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface, - FreeSurface, FFTImplicitFreeSurfaceSolver, PCGImplicitFreeSurfaceSolver, MatrixImplicitFreeSurfaceSolver, @@ -16,7 +15,6 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: using Oceananigans.Grids: with_halo - function set_simple_divergent_velocity!(model) # Create a divergent velocity grid = model.grid diff --git a/test/test_split_explicit_free_surface_solver.jl b/test/test_split_explicit_free_surface_solver.jl index b4e7a53b50..2786bf70c9 100644 --- a/test/test_split_explicit_free_surface_solver.jl +++ b/test/test_split_explicit_free_surface_solver.jl @@ -1,10 +1,15 @@ include("dependencies_for_runtests.jl") -using Oceananigans.Models.HydrostaticFreeSurfaceModels -import Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitFreeSurface -import Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitState, SplitExplicitAuxiliaryFields, SplitExplicitSettings, iterate_split_explicit! -using Oceananigans.Models.HydrostaticFreeSurfaceModels: constant_averaging_kernel -using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calculate_adaptive_settings +using Oceananigans.Models.HydrostaticFreeSurfaceModels +using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, + calculate_adaptive_settings, + constant_averaging_kernel, + materialize_free_surface, + SplitExplicitFreeSurface, + SplitExplicitState, + SplitExplicitAuxiliaryFields, + SplitExplicitSettings, + iterate_split_explicit! @testset "Split-Explicit Dynamics" begin @@ -19,10 +24,10 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc grid = RectilinearGrid(arch, FT; topology, size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (-Lz, 0), - halo=(1, 1, 1)) + halo = (1, 1, 1)) - settings = SplitExplicitSettings(eltype(grid); substeps = 200, averaging_kernel = constant_averaging_kernel) - sefs = SplitExplicitFreeSurface(grid; settings) + sefs = SplitExplicitFreeSurface(substeps = 200, averaging_kernel = constant_averaging_kernel) + sefs = materialize_free_surface(sefs, nothing, grid) sefs.η .= 0 @@ -35,9 +40,9 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc η₀(x, y, z) = sin(x) set!(η, η₀) - - Nsubsteps = calculate_substeps(settings.substepping, 1) - fractional_Δt, weights = calculate_adaptive_settings(settings.substepping, Nsubsteps) # barotropic time step in fraction of baroclinic step and averaging weights + + Nsubsteps = calculate_substeps(sefs.settings.substepping, 1) + fractional_Δt, weights = calculate_adaptive_settings(sefs.settings.substepping, Nsubsteps) # barotropic time step in fraction of baroclinic step and averaging weights iterate_split_explicit!(sefs, grid, Δτ, weights, Val(1)) @@ -59,10 +64,10 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc Nt = floor(Int, T / Δτ) Δτ_end = T - Nt * Δτ - settings = SplitExplicitSettings(; substeps = Nt, averaging_kernel = constant_averaging_kernel) - sefs = sefs(settings) + sefs = SplitExplicitFreeSurface(substeps = Nt, averaging_kernel = constant_averaging_kernel) + sefs = materialize_free_surface(sefs, nothing, grid) - # set!(η, f(x,y)) + # set!(η, f(x, y)) η₀(x, y, z) = sin(x) set!(η, η₀) U₀(x, y, z) = 0 @@ -76,13 +81,13 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc Gᵁ .= 0 Gⱽ .= 0 - weights = settings.substepping.averaging_weights - - for i in 1:Nt + weights = sefs.settings.substepping.averaging_weights + + for _ in 1:Nt iterate_split_explicit!(sefs, grid, Δτ, weights, Val(1)) end iterate_split_explicit!(sefs, grid, Δτ_end, weights, Val(1)) - + U_computed = Array(deepcopy(interior(U))) η_computed = Array(deepcopy(interior(η))) set!(η, η₀) @@ -91,13 +96,11 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc η_exact = Array(deepcopy(interior(η))) @test maximum(abs.(U_computed - U_exact)) < 1e-3 - @show maximum(abs.(η_computed)) - @show maximum(abs.(η_exact)) @test maximum(abs.(η_computed - η_exact)) < max(100eps(FT), 1e-6) end - settings = SplitExplicitSettings(eltype(grid); substeps = 200, averaging_kernel = constant_averaging_kernel) - sefs = SplitExplicitFreeSurface(grid; settings) + sefs = SplitExplicitFreeSurface(substeps = 200, averaging_kernel = constant_averaging_kernel) + sefs = materialize_free_surface(sefs, nothing, grid) sefs.η .= 0 @@ -132,7 +135,7 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc fractional_Δt, weights = calculate_adaptive_settings(settings.substepping, Nsubsteps) # barotropic time step in fraction of baroclinic step and averaging weights for step in 1:Nsubsteps - iterate_split_explicit!(sefs, grid, Δτ, weights, Val(1)) + iterate_split_explicit!(sefs, grid, Δτ, weights, Val(1)) end U_computed = Array(deepcopy(interior(U))) @@ -166,7 +169,9 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc Nt = floor(Int, T / Δτ) Δτ_end = T - Nt * Δτ - sefs = SplitExplicitFreeSurface(grid) + sefs = SplitExplicitFreeSurface(substeps=200) + sefs = materialize_free_surface(sefs, nothing, grid) + state = sefs.state auxiliary = sefs.auxiliary U, V, η̅, U̅, V̅ = state.U, state.V, state.η̅, state.U̅, state.V̅ @@ -174,7 +179,7 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc η = sefs.η g = sefs.gravitational_acceleration - # set!(η, f(x,y)) k² = ω² + # set!(η, f(x, y)) k² = ω² gu_c = 1 gv_c = 2 η₀(x, y, z) = sin(kx * x) * sin(ky * y) + 1 @@ -190,14 +195,14 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc Gᵁ .= gu_c Gⱽ .= gv_c - settings = SplitExplicitSettings(eltype(grid); substeps = Nt + 1, averaging_kernel = constant_averaging_kernel) + settings = SplitExplicitSettings(grid; substeps = Nt + 1, averaging_kernel = constant_averaging_kernel) sefs = sefs(settings) weights = settings.substepping.averaging_weights for i in 1:Nt - iterate_split_explicit!(sefs, grid, Δτ, weights, Val(1)) + iterate_split_explicit!(sefs, grid, Δτ, weights, Val(1)) end - iterate_split_explicit!(sefs, grid, Δτ_end, weights, Val(1)) + iterate_split_explicit!(sefs, grid, Δτ_end, weights, Val(1)) η_mean_after = mean(Array(interior(η))) @@ -241,4 +246,4 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: calculate_substeps, calc end end # end of architecture loop end # end of float type loop -end # end of testset loop \ No newline at end of file +end # end of testset loop diff --git a/test/test_split_explicit_vertical_integrals.jl b/test/test_split_explicit_vertical_integrals.jl index 6016053926..ca7713ad2b 100644 --- a/test/test_split_explicit_vertical_integrals.jl +++ b/test/test_split_explicit_vertical_integrals.jl @@ -1,43 +1,51 @@ include("dependencies_for_runtests.jl") + using Oceananigans.Models.HydrostaticFreeSurfaceModels using Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitFreeSurface -using Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitState, SplitExplicitAuxiliaryFields, SplitExplicitSettings -using Oceananigans.Models.HydrostaticFreeSurfaceModels: compute_barotropic_mode!, barotropic_split_explicit_corrector!, initialize_free_surface_state! +using Oceananigans.Models.HydrostaticFreeSurfaceModels: SplitExplicitState, + SplitExplicitAuxiliaryFields, + SplitExplicitSettings +using Oceananigans.Models.HydrostaticFreeSurfaceModels: compute_barotropic_mode!, + barotropic_split_explicit_corrector!, + initialize_free_surface_state! @testset "Barotropic Kernels" begin for arch in archs FT = Float64 topology = (Periodic, Periodic, Bounded) - Nx = Ny = 16 * 8 - Nz = 32 - Nx = 128 - Ny = 64 + Nx, Ny, Nz = 128, 64, 32 Lx = Ly = Lz = 2π + grid = RectilinearGrid(arch, topology = topology, size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (-Lz, 0)) - tmp = SplitExplicitFreeSurface(; substeps = 200) + tmp = SplitExplicitFreeSurface(substeps = 200) + sefs = SplitExplicitState(grid, tmp.settings.timestepper) sefs = SplitExplicitAuxiliaryFields(grid) - sefs = SplitExplicitFreeSurface(grid) + sefs = SplitExplicitFreeSurface(substeps = 200) + sefs = materialize_free_surface(sefs, nothing, grid) state = sefs.state auxiliary = sefs.auxiliary U, V, η̅, U̅, V̅ = state.U, state.V, state.η̅, state.U̅, state.V̅ Gᵁ, Gⱽ = auxiliary.Gᵁ, auxiliary.Gⱽ - u = Field{Face,Center,Center}(grid) - v = Field{Center,Face,Center}(grid) + u = Field{Face, Center, Center}(grid) + v = Field{Center, Face, Center}(grid) @testset "Average to zero" begin # set equal to something else η̅ .= U̅ .= V̅ .= 1.0 + # now set equal to zero initialize_free_surface_state!(sefs.state, sefs.η, sefs.settings.timestepper) - # don't forget the ghost points + + # don't forget the halo points fill_halo_regions!(η̅) fill_halo_regions!(U̅) fill_halo_regions!(V̅) + # check @test all(Array(η̅.data.parent) .== 0.0) @test all(Array(U̅.data.parent .== 0.0)) @@ -102,22 +110,21 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: compute_barotropic_mode! # Test 4: Test Barotropic Correction FT = Float64 topology = (Periodic, Periodic, Bounded) - Nx = Ny = 16 * 8 - Nz = 32 - Nx = 128 - Ny = 64 + Nx, Ny, Nz = 128, 64, 32 Lx = Ly = Lz = 2π + grid = RectilinearGrid(arch, topology = topology, size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (-Lz, 0)) - sefs = SplitExplicitFreeSurface(grid) + sefs = SplitExplicitFreeSurface(grid, cfl=0.7) + sefs = materialize_free_surface(sefs, nothing, grid) state = sefs.state auxiliary = sefs.auxiliary U, V, η̅, U̅, V̅ = state.U, state.V, state.η̅, state.U̅, state.V̅ Gᵁ, Gⱽ = auxiliary.Gᵁ, auxiliary.Gⱽ - u = Field{Face,Center,Center}(grid) - v = Field{Center,Face,Center}(grid) + u = Field{Face, Center, Center}(grid) + v = Field{Center, Face, Center}(grid) u_corrected = similar(u) v_corrected = similar(v) diff --git a/validation/distributed_simulations/distributed_geostrophic_adjustment.jl b/validation/distributed_simulations/distributed_geostrophic_adjustment.jl index 9659b82246..6bcbc5fd8a 100644 --- a/validation/distributed_simulations/distributed_geostrophic_adjustment.jl +++ b/validation/distributed_simulations/distributed_geostrophic_adjustment.jl @@ -43,7 +43,7 @@ grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bottom), true) model = HydrostaticFreeSurfaceModel(; grid, coriolis = FPlane(f=1e-4), - free_surface = SplitExplicitFreeSurface(; substeps=10)) + free_surface = SplitExplicitFreeSurface(grid; substeps=10)) gaussian(x, L) = exp(-x^2 / 2L^2) diff --git a/validation/distributed_simulations/distributed_hydrostatic_turbulence.jl b/validation/distributed_simulations/distributed_hydrostatic_turbulence.jl index 1358bb59cf..44f1d2b0ab 100644 --- a/validation/distributed_simulations/distributed_hydrostatic_turbulence.jl +++ b/validation/distributed_simulations/distributed_hydrostatic_turbulence.jl @@ -23,7 +23,7 @@ function run_simulation(nx, ny, arch; topology = (Periodic, Periodic, Bounded)) model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = VectorInvariant(vorticity_scheme=WENO(order=9)), - free_surface = SplitExplicitFreeSurface(substeps=10), + free_surface = SplitExplicitFreeSurface(grid, substeps=10), tracer_advection = WENO(), buoyancy = nothing, coriolis = FPlane(f = 1), diff --git a/validation/field_time_series_boundary_conditions/seasonal_baroclinic_channel.jl b/validation/field_time_series_boundary_conditions/seasonal_baroclinic_channel.jl index fb90a9307a..9a9cebdde3 100644 --- a/validation/field_time_series_boundary_conditions/seasonal_baroclinic_channel.jl +++ b/validation/field_time_series_boundary_conditions/seasonal_baroclinic_channel.jl @@ -82,7 +82,7 @@ tracer_advection = WENO() buoyancy = SeawaterBuoyancy(equation_of_state = LinearEquationOfState(), constant_salinity = 35) -free_surface = SplitExplicitFreeSurface(; grid, cfl = 0.7) +free_surface = SplitExplicitFreeSurface(grid, cfl = 0.7) coriolis = HydrostaticSphericalCoriolis() diff --git a/validation/implicit_free_surface/geostrophic_adjustment_test.jl b/validation/implicit_free_surface/geostrophic_adjustment_test.jl index bedc84db7e..45d38d8987 100644 --- a/validation/implicit_free_surface/geostrophic_adjustment_test.jl +++ b/validation/implicit_free_surface/geostrophic_adjustment_test.jl @@ -107,23 +107,12 @@ function run_and_analyze(simulation) return (ηarr, varr, uarr) end -# fft_based_free_surface = ImplicitFreeSurface() -pcg_free_surface = ImplicitFreeSurface(solver_method = :PreconditionedConjugateGradient); -matrix_free_surface = ImplicitFreeSurface(solver_method = :HeptadiagonalIterativeSolver); - -using Oceananigans.Models.HydrostaticFreeSurfaceModels: averaging_shape_function, - averaging_cosine_function, - averaging_fixed_function, +using Oceananigans.Models.HydrostaticFreeSurfaceModels: averaging_shape_function, ForwardBackwardScheme, AdamsBashforth3Scheme # @inline new_function(t) = averaging_shape_function(t; r = 0.15766, p = 2, q = 2) -@inline new_function(t) = t >= 1/2 && t <= 3/2 ? 1.0 : 0.0 - -splitexplicit_free_surface = SplitExplicitFreeSurface(substeps = 10, - averaging_weighting_function = averaging_shape_function, - timestepper = AdamsBashforth3Scheme()) -explicit_free_surface = ExplicitFreeSurface() +@inline new_function(t) = t ≥ 1/2 && t ≤ 3/2 ? 1.0 : 0.0 topology_types = [(Bounded, Periodic, Bounded), (Periodic, Periodic, Bounded)] topology_types = [topology_types[2]] @@ -131,14 +120,12 @@ topology_types = [topology_types[2]] archs = [Oceananigans.CPU()] #, Oceananigans.GPU()] archs = [archs[1]] -free_surfaces = [splitexplicit_free_surface, matrix_free_surface, explicit_free_surface] #, matrix_free_surface]; - simulations = [geostrophic_adjustment_simulation(free_surface, topology_type, multi_region, arch = arch) for (free_surface, multi_region) in zip(free_surfaces, (true, false, false)), topology_type in topology_types, - arch in archs]; + arch in archs] -data = [run_and_analyze(sim) for sim in simulations]; +data = [run_and_analyze(sim) for sim in simulations] using GLMakie using JLD2 @@ -153,6 +140,17 @@ grid = RectilinearGrid(arch, x = (0, Lh), y = (0, Lh), z = (-Lz, 0), topology = topology) +explicit_free_surface = ExplicitFreeSurface() +fft_based_free_surface = ImplicitFreeSurface(solver_method = :FastFourierTransform) +pcg_free_surface = ImplicitFreeSurface(solver_method = :PreconditionedConjugateGradient) +matrix_free_surface = ImplicitFreeSurface(solver_method = :HeptadiagonalIterativeSolver) +splitexplicit_free_surface = SplitExplicitFreeSurface(grid, + substeps = 10, + averaging_kernel = averaging_shape_function, + timestepper = AdamsBashforth3Scheme()) + +free_surfaces = [splitexplicit_free_surface, matrix_free_surface, explicit_free_surface] + x = grid.xᶜᵃᵃ[1:grid.Nx] xf = grid.xᶠᵃᵃ[1:grid.Nx+1] y = grid.yᵃᶜᵃ[1:grid.Ny] @@ -220,7 +218,3 @@ GLMakie.record(fig, "free_surface_bounded.mp4", 1:1000, framerate = 12) do i @info "Plotting iteration $i of 600..." iter[] = i end - - - -