diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl index 588235882e..7bd3ee6337 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl @@ -201,7 +201,7 @@ function HydrostaticFreeSurfaceModel(; grid, free_surface, forcing, closure, particles, biogeochemistry, velocities, tracers, pressure, diffusivity_fields, timestepper, auxiliary_fields) - update_state!(model) + update_state!(model; compute_tendencies = false) return model end diff --git a/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl b/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl index 0238f876cb..f216d312c7 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/set_hydrostatic_free_surface_model.jl @@ -61,7 +61,7 @@ model.velocities.u end initialize!(model) - update_state!(model) + update_state!(model; compute_tendencies = false) return nothing end diff --git a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl index 8429887ea7..d4d3dec8d8 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl @@ -61,7 +61,7 @@ compute_free_surface_tendency!(::SingleColumnGrid, ::SplitExplicitFreeSurfaceHFS # Fast state update and halo filling -function update_state!(model::HydrostaticFreeSurfaceModel, grid::SingleColumnGrid, callbacks) +function update_state!(model::HydrostaticFreeSurfaceModel, grid::SingleColumnGrid, callbacks; compute_tendencies = true) fill_halo_regions!(prognostic_fields(model), model.clock, fields(model)) @@ -78,6 +78,9 @@ function update_state!(model::HydrostaticFreeSurfaceModel, grid::SingleColumnGri end update_biogeochemical_state!(model.biogeochemistry, model) + + compute_tendencies && + @apply_regionally compute_tendencies!(model, callbacks) return nothing end 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 fe94fbe3fa..cf6b20eade 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/update_hydrostatic_free_surface_model_state.jl @@ -19,11 +19,12 @@ compute_auxiliary_fields!(auxiliary_fields) = Tuple(compute!(a) for a in auxilia # single column models. """ - update_state!(model::HydrostaticFreeSurfaceModel, callbacks=[]) + update_state!(model::HydrostaticFreeSurfaceModel, callbacks=[]; compute_tendencies = true) Update peripheral aspects of the model (auxiliary fields, halo regions, diffusivities, hydrostatic pressure) to the current model state. If `callbacks` are provided (in an array), -they are called in the end. +they are called in the end. Finally, the tendencies for the new time-step are computed if +`compute_tendencies = true`. """ update_state!(model::HydrostaticFreeSurfaceModel, callbacks=[]; compute_tendencies = true) = update_state!(model, model.grid, callbacks; compute_tendencies) diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl index 688f834f3a..0b42767b9a 100644 --- a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl +++ b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl @@ -230,8 +230,8 @@ function NonhydrostaticModel(; grid, pressures, diffusivity_fields, timestepper, pressure_solver, immersed_boundary, auxiliary_fields) - update_state!(model) - + update_state!(model; compute_tendencies = false) + return model end diff --git a/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl index e8011c9058..83dcda0807 100644 --- a/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl +++ b/src/Models/NonhydrostaticModels/set_nonhydrostatic_model.jl @@ -47,13 +47,13 @@ function set!(model::NonhydrostaticModel; enforce_incompressibility=true, kwargs # Apply a mask foreach(mask_immersed_field!, model.tracers) foreach(mask_immersed_field!, model.velocities) - update_state!(model) + update_state!(model; compute_tendencies = false) if enforce_incompressibility FT = eltype(model.grid) calculate_pressure_correction!(model, one(FT)) pressure_correct_velocities!(model, one(FT)) - update_state!(model) + update_state!(model; compute_tendencies = false) end return nothing diff --git a/src/Models/ShallowWaterModels/set_shallow_water_model.jl b/src/Models/ShallowWaterModels/set_shallow_water_model.jl index b68f279a04..b836106d12 100644 --- a/src/Models/ShallowWaterModels/set_shallow_water_model.jl +++ b/src/Models/ShallowWaterModels/set_shallow_water_model.jl @@ -14,7 +14,7 @@ function set!(model::ShallowWaterModel; kwargs...) set!(ϕ, value) end - update_state!(model) + update_state!(model; compute_tendencies = false) return nothing end diff --git a/src/Models/ShallowWaterModels/shallow_water_model.jl b/src/Models/ShallowWaterModels/shallow_water_model.jl index 9c5778ac14..7e51210703 100644 --- a/src/Models/ShallowWaterModels/shallow_water_model.jl +++ b/src/Models/ShallowWaterModels/shallow_water_model.jl @@ -205,7 +205,7 @@ function ShallowWaterModel(; timestepper, formulation) - update_state!(model) + update_state!(model; compute_tendencies = false) return model end diff --git a/src/TimeSteppers/quasi_adams_bashforth_2.jl b/src/TimeSteppers/quasi_adams_bashforth_2.jl index 891d8fbb99..ec047f4a61 100644 --- a/src/TimeSteppers/quasi_adams_bashforth_2.jl +++ b/src/TimeSteppers/quasi_adams_bashforth_2.jl @@ -57,7 +57,7 @@ reset!(timestepper::QuasiAdamsBashforth2TimeStepper) = nothing ##### """ - time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false, compute_tendencies=true) + time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; euler=false) Step forward `model` one time step `Δt` with a 2nd-order Adams-Bashforth method and pressure-correction substep. Setting `euler=true` will take a forward Euler time step. @@ -71,14 +71,10 @@ The steps of the Quasi-Adams-Bashforth second-order (AB2) algorithm are: 4. Correct the velocities based on the results of step 3. 5. Store the old tendencies. 6. Update the model state. -7. If `compute_tendencies == true`, compute the tendencies for the next time step. - -!!! warning "`compute_tendencies` kwarg" - If `compute_tendencies == false` then the new tendencies at the last step are _not_ calculated! - Setting `compute_tendencies == false` is not recommended _except_ for debugging purposes. +7. Compute tendencies for the next time step """ function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt; - callbacks=[], euler=false, compute_tendencies=true) + callbacks=[], euler=false) Δt == 0 && @warn "Δt == 0 may cause model blowup!" @@ -113,7 +109,7 @@ function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt end # Be paranoid and update state at iteration 0 - model.clock.iteration == 0 && update_state!(model, callbacks) + model.clock.iteration == 0 && update_state!(model, callbacks; compute_tendencies=true) ab2_step!(model, Δt) # full step for tracers, fractional step for velocities. @@ -124,7 +120,7 @@ function time_step!(model::AbstractModel{<:QuasiAdamsBashforth2TimeStepper}, Δt calculate_pressure_correction!(model, Δt) @apply_regionally correct_velocities_and_store_tendencies!(model, Δt) - update_state!(model, callbacks; compute_tendencies) + update_state!(model, callbacks; compute_tendencies=true) step_lagrangian_particles!(model, Δt) # Return χ to initial value diff --git a/src/TimeSteppers/runge_kutta_3.jl b/src/TimeSteppers/runge_kutta_3.jl index c54023e65a..0bc5c34608 100644 --- a/src/TimeSteppers/runge_kutta_3.jl +++ b/src/TimeSteppers/runge_kutta_3.jl @@ -78,11 +78,11 @@ The 3rd-order Runge-Kutta method takes three intermediate substep stages to achieve a single timestep. A pressure correction step is applied at each intermediate stage. """ -function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbacks=[], compute_tendencies = true) +function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbacks=[]) Δt == 0 && @warn "Δt == 0 may cause model blowup!" # Be paranoid and update state at iteration 0, in case run! is not used: - model.clock.iteration == 0 && update_state!(model, callbacks) + model.clock.iteration == 0 && update_state!(model, callbacks; compute_tendencies = true) γ¹ = model.timestepper.γ¹ γ² = model.timestepper.γ² @@ -111,7 +111,7 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac pressure_correct_velocities!(model, first_stage_Δt) store_tendencies!(model) - update_state!(model, callbacks) + update_state!(model, callbacks; compute_tendencies = true) step_lagrangian_particles!(model, first_stage_Δt) # @@ -127,7 +127,7 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac pressure_correct_velocities!(model, second_stage_Δt) store_tendencies!(model) - update_state!(model, callbacks) + update_state!(model, callbacks; compute_tendencies = true) step_lagrangian_particles!(model, second_stage_Δt) # @@ -148,7 +148,7 @@ function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; callbac calculate_pressure_correction!(model, third_stage_Δt) pressure_correct_velocities!(model, third_stage_Δt) - update_state!(model, callbacks; compute_tendencies) + update_state!(model, callbacks; compute_tendencies = true) step_lagrangian_particles!(model, third_stage_Δt) return nothing diff --git a/test/test_ensemble_hydrostatic_free_surface_models.jl b/test/test_ensemble_hydrostatic_free_surface_models.jl index 17736968be..95ab1b701c 100644 --- a/test/test_ensemble_hydrostatic_free_surface_models.jl +++ b/test/test_ensemble_hydrostatic_free_surface_models.jl @@ -4,6 +4,41 @@ using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize, Slic using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity const CAVD = ConvectiveAdjustmentVerticalDiffusivity +@testset "`HydrostaticFreeSurfaceModel` using a `SingleColumnGrid`" begin + + Nz = 3 + Hz = 1 + single_column_topology = (Flat, Flat, Bounded) + periodic_topology = (Periodic, Periodic, Bounded) + + single_column_grid = RectilinearGrid(; size=Nz, z=(-1, 0), topology = single_column_topology, halo=Hz) + periodic_grid = RectilinearGrid(; size=(1, 1, Nz), x = (0, 1), y = (0, 1), z=(-1, 0), topology = periodic_topology, halo=(1, 1, Hz)) + coriolis = FPlane(f=0.2) + closure = CAVD(background_κz=1.0) + + Δt = 0.01 + + model_kwargs = (; tracers=:c, buoyancy=nothing, closure, coriolis) + simulation_kwargs = (; Δt, stop_iteration=100) + + sic_model = HydrostaticFreeSurfaceModel(; grid = single_column_grid, model_kwargs...) + per_model = HydrostaticFreeSurfaceModel(; grid = periodic_grid, model_kwargs...) + + set!(sic_model, c = z -> exp(-z^2), u = 1, v = 1) + set!(per_model, c = (x, y, z) -> exp(-z^2), u = 1, v = 1) + + sic_simulation = Simulation(sic_model; simulation_kwargs...) + per_simulation = Simulation(per_model; simulation_kwargs...) + run!(sic_simulation) + run!(per_simulation) + + @info "Testing Single column grid results..." + + @test all(sic_model.velocities.u.data[1, 1, :] .≈ per_model.velocities.u.data[1, 1, :]) + @test all(sic_model.velocities.v.data[1, 1, :] .≈ per_model.velocities.v.data[1, 1, :]) + @test all(sic_model.tracers.c.data[1, 1, :] .≈ per_model.tracers.c.data[1, 1, :]) +end + @testset "Ensembles of `HydrostaticFreeSurfaceModel` with different closures" begin Nz = 16 diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index 6f5334ec87..50fd4ff90f 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -99,11 +99,12 @@ function run_first_AB2_time_step_tests(arch, FT) buoyancy = SeawaterBuoyancy(), tracers = (:T, :S)) - # Test that GT = 1 after model construction (note: this computes tendencies) + # Test that GT = 0 after model construction + # (note: model construction does not computes tendencies) @test all(interior(model.timestepper.Gⁿ.u) .≈ 0) @test all(interior(model.timestepper.Gⁿ.v) .≈ 0) @test all(interior(model.timestepper.Gⁿ.w) .≈ 0) - @test all(interior(model.timestepper.Gⁿ.T) .≈ 1) + @test all(interior(model.timestepper.Gⁿ.T) .≈ 0) @test all(interior(model.timestepper.Gⁿ.S) .≈ 0) # Test that T = 1 after 1 time step and that AB2 actually reduced to forward Euler.