From 75bb068211c51671a0ea721e333c13870ea18038 Mon Sep 17 00:00:00 2001 From: Brandon Allen Date: Wed, 22 Jul 2020 15:05:08 -0400 Subject: [PATCH] Add traits to modify HBModel and SWModel for coupling drive shallow water model by HB model add coupling traits, use to replace forcing traits in shallow water model move coupling for SWModel to file in SplitExplicit add coupling to HBModel, get fully coupled model to converge to analytic solution Add hydrostatic spindown tests with actual multi-rate happening remove old integral models get flow deviation working with new reshape&broadcast setup (on cpu). Remove stub functions for partial coupling use Coupling trait for dispatch instead of boolean in SplitExplicitMethod. use abbreviated model names Update to new functions from #1280 Co-Authored-By: Jean-Michel Campin --- .../ODESolvers/SplitExplicitMethod.jl | 113 ++++++++---------- .../HydrostaticBoussinesqModel.jl | 108 +++++++++++------ src/Ocean/Ocean.jl | 8 ++ src/Ocean/ShallowWater/ShallowWaterModel.jl | 30 +++-- .../HydrostaticBoussinesqCoupling.jl | 81 +++++++++++++ .../SplitExplicit/ShallowWaterCoupling.jl | 7 ++ src/Ocean/SplitExplicit/SplitExplicitModel.jl | 34 ++++++ .../ShallowWater/2D_hydrostatic_spindown.jl | 2 + test/Ocean/ShallowWater/GyreDriver.jl | 9 +- test/Ocean/ShallowWater/GyreInABox.jl | 1 + .../SplitExplicit/hydrostatic_spindown.jl | 67 +++++++++-- .../test_vertical_integral_model.jl | 2 + 12 files changed, 344 insertions(+), 118 deletions(-) create mode 100644 src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl create mode 100644 src/Ocean/SplitExplicit/ShallowWaterCoupling.jl diff --git a/src/Numerics/ODESolvers/SplitExplicitMethod.jl b/src/Numerics/ODESolvers/SplitExplicitMethod.jl index 13b0e1bd455..fa369e3bb8c 100644 --- a/src/Numerics/ODESolvers/SplitExplicitMethod.jl +++ b/src/Numerics/ODESolvers/SplitExplicitMethod.jl @@ -36,8 +36,6 @@ mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver dt::RT "time" t::RT - "flag to couple the models" - coupled "storage for transfer tendency" dQ2fast::MSA @@ -47,7 +45,6 @@ mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver Q = nothing; dt = getdt(slow_solver), t0 = slow_solver.t, - coupled = true, ) where {AT <: AbstractArray} SS = typeof(slow_solver) FS = typeof(fast_solver) @@ -62,7 +59,6 @@ mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver fast_solver, RT(dt), RT(t0), - coupled, dQ2fast, ) end @@ -97,53 +93,28 @@ function dostep!( # Initialize fast model and tendency adjustment # before evalution of slow mode - if split.coupled - initialize_states!( - slow_bl, - fast_bl, - slow.rhs!, - fast.rhs!, - Qslow, - Qfast, - ) - - # Evaluate the slow mode - # --> save tendency for the fast - slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false) - - # vertically integrate slow tendency to advance fast equation - # and use vertical mean for slow model (negative source) - # ---> work with dQ2fast as input - tendency_from_slow_to_fast!( - slow_bl, - fast_bl, - slow.rhs!, - fast.rhs!, - Qslow, - Qfast, - dQ2fast, - ) - end + initialize_states!(slow_bl, fast_bl, slow.rhs!, fast.rhs!, Qslow, Qfast) + + # Evaluate the slow mode + # --> save tendency for the fast + slow.rhs!(dQ2fast, Qslow, param, slow_stage_time, increment = false) + + # vertically integrate slow tendency to advance fast equation + # and use vertical mean for slow model (negative source) + # ---> work with dQ2fast as input + tendency_from_slow_to_fast!( + slow_bl, + fast_bl, + slow.rhs!, + fast.rhs!, + Qslow, + Qfast, + dQ2fast, + ) # Compute (and RK update) slow tendency slow.rhs!(dQslow, Qslow, param, slow_stage_time, increment = true) - # Update (RK-stage) slow state - event = Event(array_device(Qslow)) - event = update!(array_device(Qslow), groupsize)( - realview(dQslow), - realview(Qslow), - slow.RKA[slow_s % length(slow.RKA) + 1], - slow.RKB[slow_s], - slow_dt, - nothing, - nothing, - nothing; - ndrange = length(realview(Qslow)), - dependencies = (event,), - ) - wait(array_device(Qslow), event) - # Fractional time for slow stage if slow_s == length(slow.RKA) γ = 1 - slow.RKC[slow_s] @@ -161,30 +132,42 @@ function dostep!( for substep in 1:nsubsteps fast_time = slow_stage_time + (substep - 1) * fast_dt dostep!(Qfast, fast, param, fast_time) - if split.coupled - cummulate_fast_solution!( - slow_bl, - fast_bl, - fast.rhs!, - Qfast, - fast_time, - fast_dt, - substep, - ) - end - end - - # reconcile slow equation using fast equation - if split.coupled - reconcile_from_fast_to_slow!( + cummulate_fast_solution!( slow_bl, fast_bl, - slow.rhs!, fast.rhs!, - Qslow, Qfast, + fast_time, + fast_dt, + substep, ) end + + # Update (RK-stage) slow state + event = Event(array_device(Qslow)) + event = update!(array_device(Qslow), groupsize)( + realview(dQslow), + realview(Qslow), + slow.RKA[slow_s % length(slow.RKA) + 1], + slow.RKB[slow_s], + slow_dt, + nothing, + nothing, + nothing; + ndrange = length(realview(Qslow)), + dependencies = (event,), + ) + wait(array_device(Qslow), event) + + # reconcile slow equation using fast equation + reconcile_from_fast_to_slow!( + slow_bl, + fast_bl, + slow.rhs!, + fast.rhs!, + Qslow, + Qfast, + ) end updatedt!(fast, fast_dt_in) diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl index 74adc59dcc9..3c32f865156 100644 --- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl +++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl @@ -6,6 +6,7 @@ using StaticArrays using LinearAlgebra: dot, Diagonal using CLIMAParameters.Planet: grav +using ..Ocean using ...VariableTemplates using ...MPIStateArrays using ...Mesh.Filters: apply! @@ -13,11 +14,11 @@ using ...Mesh.Grids: VerticalDirection using ...Mesh.Geometry using ...DGMethods using ...DGMethods.NumericalFluxes -import ...DGMethods.NumericalFluxes: update_penalty! using ...DGMethods.NumericalFluxes: RusanovNumericalFlux using ...BalanceLaws -using ...BalanceLaws: number_states +import ..Ocean: coriolis_parameter! +import ...DGMethods.NumericalFluxes: update_penalty! import ...BalanceLaws: vars_state, init_state_prognostic!, @@ -67,9 +68,10 @@ fₒ = first coriolis parameter (constant term) HydrostaticBoussinesqModel(problem) """ -struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw +struct HydrostaticBoussinesqModel{C, PS, P, T} <: BalanceLaw param_set::PS problem::P + coupling::C ρₒ::T cʰ::T cᶻ::T @@ -83,7 +85,8 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw β::T function HydrostaticBoussinesqModel{FT}( param_set::PS, - problem; + problem::P; + coupling::C = Uncoupled(), ρₒ = FT(1000), # kg / m^3 cʰ = FT(0), # m/s cᶻ = FT(0), # m/s @@ -95,10 +98,11 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw κᶜ = FT(1e-1), # m^2 / s # diffusivity for convective adjustment fₒ = FT(1e-4), # Hz β = FT(1e-11), # Hz / m - ) where {FT <: AbstractFloat, PS} - return new{PS, typeof(problem), FT}( + ) where {FT <: AbstractFloat, PS, P, C} + return new{C, PS, P, FT}( param_set, problem, + coupling, ρₒ, cʰ, cᶻ, @@ -189,6 +193,7 @@ these are just copies in our model function vars_state(m::HBModel, ::Gradient, T) @vars begin ∇u::SVector{2, T} + ∇uᵈ::SVector{2, T} ∇θ::T end end @@ -207,9 +212,23 @@ this computation is done pointwise at each nodal point - `t`: time, not used """ @inline function compute_gradient_argument!(m::HBModel, G::Vars, Q::Vars, A, t) - G.∇u = Q.u G.∇θ = Q.θ + velocity_gradient_argument!(m, m.coupling, G, Q, A, t) + + return nothing +end + +@inline function velocity_gradient_argument!( + m::HBModel, + ::Uncoupled, + G, + Q, + A, + t, +) + G.∇u = Q.u + return nothing end @@ -252,8 +271,7 @@ this computation is done pointwise at each nodal point # store ∇ʰu for continuity equation (convert gradient to divergence) D.∇ʰu = G.∇u[1, 1] + G.∇u[2, 2] - ν = viscosity_tensor(m) - D.ν∇u = -ν * G.∇u + velocity_gradient_flux!(m, m.coupling, D, G, Q, A, t) κ = diffusivity_tensor(m, G.∇θ[3]) D.κ∇θ = -κ * G.∇θ @@ -261,6 +279,13 @@ this computation is done pointwise at each nodal point return nothing end +@inline function velocity_gradient_flux!(m::HBModel, ::Uncoupled, D, G, Q, A, t) + ν = viscosity_tensor(m) + D.ν∇u = -ν * G.∇u + + return nothing +end + """ viscosity_tensor(::HBModel) @@ -422,38 +447,44 @@ t -> time, not used t::Real, direction, ) - FT = eltype(Q) - _grav::FT = grav(m.param_set) @inbounds begin - u = Q.u # Horizontal components of velocity - η = Q.η - θ = Q.θ - w = A.w # vertical velocity - pkin = A.pkin + # ∇h • (g η) + hydrostatic_pressure!(m, m.coupling, F, Q, A, t) - v = @SVector [u[1], u[2], w] + # ∇h • (- ∫(αᵀ θ)) + pkin = A.pkin Iʰ = @SMatrix [ 1 -0 -0 1 -0 -0 ] - - # ∇h • (g η) - F.u += _grav * η * Iʰ - - # ∇h • (- ∫(αᵀ θ)) - F.u += _grav * pkin * Iʰ + F.u += grav(m.param_set) * pkin * Iʰ # ∇h • (v ⊗ u) # F.u += v * u' # ∇ • (u θ) + θ = Q.θ + v = @SVector [Q.u[1], Q.u[2], A.w] F.θ += v * θ end return nothing end +@inline function hydrostatic_pressure!(m::HBModel, ::Uncoupled, F, Q, A, t) + η = Q.η + Iʰ = @SMatrix [ + 1 -0 + -0 1 + -0 -0 + ] + + F.u += grav(m.param_set) * η * Iʰ + + return nothing +end + """ flux_second_order!(::HBModel) @@ -512,22 +543,25 @@ end t::Real, direction, ) - @inbounds begin - u, v = Q.u # Horizontal components of velocity - wz0 = A.wz0 + wz0 = A.wz0 + S.η += wz0 - # f × u - f = coriolis_force(m, A.y) - S.u -= @SVector [-f * v, f * u] + coriolis_force!(m, m.coupling, S, Q, A, t) - S.η += wz0 - end + return nothing +end + +@inline function coriolis_force!(m::HBModel, ::Uncoupled, S, Q, A, t) + # f × u + f = coriolis_parameter(m, A.y) + u, v = Q.u # Horizontal components of velocity + S.u -= @SVector [-f * v, f * u] return nothing end """ - coriolis_force(::HBModel) + coriolis_parameter(::HBModel) northern hemisphere coriolis @@ -535,7 +569,7 @@ northern hemisphere coriolis - `m`: model object to dispatch on and get coriolis parameters - `y`: y-coordinate in the box """ -@inline coriolis_force(m::HBModel, y) = m.fₒ + m.β * y +@inline coriolis_parameter(m::HBModel, y) = m.fₒ + m.β * y """ wavespeed(::HBModel) @@ -595,9 +629,13 @@ function update_auxiliary_state!( apply!(Q, (:θ,), dg.grid, exp_filter, direction = VerticalDirection()) end + compute_flow_deviation!(dg, m, m.coupling, Q, t) + return true end +@inline compute_flow_deviation!(dg, ::HBModel, ::Uncoupled, _...) = nothing + """ update_auxiliary_state_gradient!(::HBModel) @@ -636,13 +674,13 @@ function update_auxiliary_state_gradient!( # We are unable to use vars (ie A.w) for this because this operation will # return a SubArray, and adapt (used for broadcasting along reshaped arrays) # has a limited recursion depth for the types allowed. - number_auxiliary = number_states(m, Auxiliary(), FT) + number_aux = number_states(m, Auxiliary(), FT) index_w = varsindex(vars_state(m, Auxiliary(), FT), :w) index_wz0 = varsindex(vars_state(m, Auxiliary(), FT), :wz0) Nq, Nqk, _, _, nelemv, nelemh, nhorzrealelem, _ = basic_grid_info(dg) # project w(z=0) down the stack - data = reshape(A.data, Nq^2, Nqk, number_auxiliary, nelemv, nelemh) + data = reshape(A.data, Nq^2, Nqk, number_aux, nelemv, nelemh) flat_wz0 = @view data[:, end:end, index_w, end:end, 1:nhorzrealelem] boxy_wz0 = @view data[:, :, index_wz0, :, 1:nhorzrealelem] boxy_wz0 .= flat_wz0 diff --git a/src/Ocean/Ocean.jl b/src/Ocean/Ocean.jl index 9dba464b5f5..ec0fcf7c052 100644 --- a/src/Ocean/Ocean.jl +++ b/src/Ocean/Ocean.jl @@ -1,5 +1,13 @@ module Ocean +export AbstractOceanCoupling, Uncoupled, Coupled + +abstract type AbstractOceanCoupling end +struct Uncoupled <: AbstractOceanCoupling end +struct Coupled <: AbstractOceanCoupling end + +function coriolis_parameter! end + include("HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl") include("ShallowWater/ShallowWaterModel.jl") include("SplitExplicit/SplitExplicitModel.jl") diff --git a/src/Ocean/ShallowWater/ShallowWaterModel.jl b/src/Ocean/ShallowWater/ShallowWaterModel.jl index f44be0ac874..5ffc4fb54ad 100644 --- a/src/Ocean/ShallowWater/ShallowWaterModel.jl +++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl @@ -6,6 +6,7 @@ using StaticArrays using LinearAlgebra: dot, Diagonal using CLIMAParameters.Planet: grav +using ..Ocean using ...VariableTemplates using ...Mesh.Geometry using ...DGMethods @@ -42,9 +43,10 @@ end abstract type AdvectionTerm end struct NonLinearAdvection <: AdvectionTerm end -struct ShallowWaterModel{PS, P, T, A, S} <: BalanceLaw +struct ShallowWaterModel{C, PS, P, T, A, S} <: BalanceLaw param_set::PS problem::P + coupling::C turbulence::T advection::A c::S @@ -156,12 +158,13 @@ function compute_gradient_flux!( α::Vars, t::Real, ) - compute_gradient_flux!(m.turbulence, σ, δ, q, α, t) + compute_gradient_flux!(m, m.turbulence, σ, δ, q, α, t) end -compute_gradient_flux!(::LinearDrag, _...) = nothing +compute_gradient_flux!(::SWModel, ::LinearDrag, _...) = nothing @inline function compute_gradient_flux!( + ::SWModel, T::ConstantViscosity, σ::Vars, δ::Grad, @@ -186,12 +189,13 @@ function flux_second_order!( α::Vars, t::Real, ) - flux_second_order!(m.turbulence, G, q, σ, α, t) + flux_second_order!(m, m.turbulence, G, q, σ, α, t) end -flux_second_order!(::LinearDrag, _...) = nothing +flux_second_order!(::SWModel, ::LinearDrag, _...) = nothing @inline function flux_second_order!( + ::SWModel, ::ConstantViscosity, G::Grad, q::Vars, @@ -210,22 +214,28 @@ end m::SWModel{P}, S::Vars, q::Vars, - diffusive::Vars, + d::Vars, α::Vars, t::Real, direction, ) where {P} - S.U += α.τ - - f = α.f - U, V = q.U + # f × u + f = A.f + U, V = Q.U S.U -= @SVector [-f * V, f * U] + forcing_term!(m, m.coupling, S, q, α, t) linear_drag!(m.turbulence, S, q, α, t) return nothing end +@inline function forcing_term!(::SWModel, ::Uncoupled, S, Q, A, t) + S.U += A.τ + + return nothing +end + linear_drag!(::ConstantViscosity, _...) = nothing @inline function linear_drag!(T::LinearDrag, S::Vars, q::Vars, α::Vars, t::Real) diff --git a/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl b/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl new file mode 100644 index 00000000000..f1c22d11d38 --- /dev/null +++ b/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl @@ -0,0 +1,81 @@ +using ..HydrostaticBoussinesq +using ...DGMethods + +import ...BalanceLaws +import ..HydrostaticBoussinesq: + viscosity_tensor, + coriolis_parameter, + velocity_gradient_argument!, + velocity_gradient_flux!, + hydrostatic_pressure!, + compute_flow_deviation! + +@inline function velocity_gradient_argument!(m::HBModel, ::Coupled, G, Q, A, t) + G.∇u = Q.u + G.∇uᵈ = A.uᵈ + + return nothing +end + +@inline function velocity_gradient_flux!(m::HBModel, ::Coupled, D, G, Q, A, t) + ν = viscosity_tensor(m) + ∇u = @SMatrix [ + G.∇uᵈ[1, 1] G.∇uᵈ[1, 2] + G.∇uᵈ[2, 1] G.∇uᵈ[2, 2] + G.∇u[3, 1] G.∇u[3, 2] + ] + D.ν∇u = -ν * ∇u + + return nothing +end + +@inline hydrostatic_pressure!(::HBModel, ::Coupled, _...) = nothing + +@inline function coriolis_force!(m::HBModel, ::Coupled, S, Q, A, t) + # f × u + f = coriolis_parameter(m, A.y) + uᵈ, vᵈ = A.uᵈ # Horizontal components of velocity + S.u -= @SVector [-f * vᵈ, f * uᵈ] + + return nothing +end + +# Compute Horizontal Flow deviation from vertical mean +@inline function compute_flow_deviation!(dg, m::HBModel, ::Coupled, Q, t) + FT = eltype(Q) + Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(dg) + + #### integrate the tendency + model_int = dg.modeldata.integral_model + integral = model_int.balance_law + update_auxiliary_state!(model_int, integral, Q, 0) + + ### properly shape MPIStateArrays + num_int = number_states(integral, ::Auxiliary, FT) + data_int = model_int.state_auxiliary.data + data_int = reshape(data_int, Nq^2, Nqk, num_int, nelemv, nelemh) + + num_aux = number_states(m, ::Auxiliary, FT) + data_aux = dg.state_auxiliary.data + data_aux = reshape(data_aux, Nq^2, Nqk, num_aux, nelemv, nelemh) + + num_state = number_states(m, ::Prognostic, FT) + data_state = reshape(Q.data, Nq^2, Nqk, num_state, nelemv, nelemh) + + ### get vars indices + index_∫u = varsindex(vars_state(integral, ::Auxiliary, FT), :(∫x)) + index_uᵈ = varsindex(vars_state(m, ::Auxiliary, FT), :uᵈ) + index_u = varsindex(vars_state(m, ::Prognostic, FT), :u) + + ### get top value (=integral over full depth) + ∫u = @view data_int[:, end:end, index_∫u, end:end, 1:nrealelemh] + uᵈ = @view data_aux[:, :, index_uᵈ, :, 1:nrealelemh] + u = @view data_state[:, :, index_u, :, 1:nrealelemh] + + ## make a copy of horizontal velocity + ## and remove vertical mean velocity + uᵈ .= u + uᵈ .-= ∫u / m.problem.H + + return nothing +end diff --git a/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl b/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl new file mode 100644 index 00000000000..9c19e161bf7 --- /dev/null +++ b/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl @@ -0,0 +1,7 @@ +import ..ShallowWater: forcing_term! + +@inline function forcing_term!(::SWModel, ::Coupled, S, Q, A, t) + S.U += A.Gᵁ + + return nothing +end diff --git a/src/Ocean/SplitExplicit/SplitExplicitModel.jl b/src/Ocean/SplitExplicit/SplitExplicitModel.jl index 1abb817ff3e..46c73fb0f83 100644 --- a/src/Ocean/SplitExplicit/SplitExplicitModel.jl +++ b/src/Ocean/SplitExplicit/SplitExplicitModel.jl @@ -2,6 +2,7 @@ module SplitExplicit using StaticArrays +using ..Ocean using ..HydrostaticBoussinesq using ..ShallowWater @@ -17,8 +18,41 @@ import ...BalanceLaws: cummulate_fast_solution!, reconcile_from_fast_to_slow! +HBModel = HydrostaticBoussinesqModel +SWModel = ShallowWaterModel + +function initialize_states!( + ::HBModel{C}, + ::SWModel{C}, + _..., +) where {C <: Uncoupled} + return nothing +end +function tendency_from_slow_to_fast!( + ::HBModel{C}, + ::SWModel{C}, + _..., +) where {C <: Uncoupled} + return nothing +end +function cummulate_fast_solution!( + ::HBModel{C}, + ::SWModel{C}, + _..., +) where {C <: Uncoupled} + return nothing +end +function reconcile_from_fast_to_slow!( + ::HBModel{C}, + ::SWModel{C}, + _..., +) where {C <: Uncoupled} + return nothing +end include("VerticalIntegralModel.jl") include("Communication.jl") +include("ShallowWaterCoupling.jl") +include("HydrostaticBoussinesqCoupling.jl") end diff --git a/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl b/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl index 5b1416b3dd5..a2541339917 100644 --- a/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl +++ b/test/Ocean/ShallowWater/2D_hydrostatic_spindown.jl @@ -7,6 +7,7 @@ using ClimateMachine.ODESolvers using ClimateMachine.Mesh.Filters using ClimateMachine.VariableTemplates using ClimateMachine.Mesh.Grids: polynomialorder +using ClimateMachine.Ocean using ClimateMachine.Ocean.ShallowWater using ClimateMachine.Mesh.Topologies @@ -102,6 +103,7 @@ function run_hydrostatic_spindown(; refDat = ()) model_2D = ShallowWaterModel( param_set, prob_2D, + Uncoupled(), ShallowWater.ConstantViscosity{FT}(5e3), nothing, FT(1), diff --git a/test/Ocean/ShallowWater/GyreDriver.jl b/test/Ocean/ShallowWater/GyreDriver.jl index 5779b8628f1..f736854c3ef 100644 --- a/test/Ocean/ShallowWater/GyreDriver.jl +++ b/test/Ocean/ShallowWater/GyreDriver.jl @@ -58,7 +58,14 @@ function setup_model(FT, stommel, linear, τₒ, fₒ, β, γ, ν, Lˣ, Lʸ, H) advection = NonLinearAdvection() end - model = ShallowWaterModel(param_set, problem, turbulence, advection, c) + model = ShallowWaterModel( + param_set, + problem, + Uncoupled(), + turbulence, + advection, + c, + ) end function shallow_init_state!( diff --git a/test/Ocean/ShallowWater/GyreInABox.jl b/test/Ocean/ShallowWater/GyreInABox.jl index bf3d81d0139..90d19f0d39d 100644 --- a/test/Ocean/ShallowWater/GyreInABox.jl +++ b/test/Ocean/ShallowWater/GyreInABox.jl @@ -8,6 +8,7 @@ using ClimateMachine.MPIStateArrays using ClimateMachine.ODESolvers using ClimateMachine.GenericCallbacks using ClimateMachine.VariableTemplates: flattenednames +using ClimateMachine.Ocean using ClimateMachine.Ocean.ShallowWater using LinearAlgebra using StaticArrays diff --git a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl index a9824031672..732781b290f 100644 --- a/test/Ocean/SplitExplicit/hydrostatic_spindown.jl +++ b/test/Ocean/SplitExplicit/hydrostatic_spindown.jl @@ -7,6 +7,7 @@ using ClimateMachine.ODESolvers using ClimateMachine.Mesh.Filters using ClimateMachine.VariableTemplates using ClimateMachine.Mesh.Grids: polynomialorder +using ClimateMachine.Ocean using ClimateMachine.Ocean.HydrostaticBoussinesq using ClimateMachine.Ocean.ShallowWater using ClimateMachine.Ocean.SplitExplicit: VerticalIntegralModel @@ -85,7 +86,11 @@ function shallow_init_aux!(m::ShallowWaterModel, p::GyreInABox, A, geom) return nothing end -function run_hydrostatic_spindown(; coupled = true, refDat = ()) +function run_hydrostatic_spindown(; + coupling = Uncoupled(), + dt_slow = 300, + refDat = (), +) mpicomm = MPI.COMM_WORLD ArrayType = ClimateMachine.array_type() @@ -130,6 +135,7 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ()) model_3D = HydrostaticBoussinesqModel{FT}( param_set, prob_3D; + coupling = coupling, cʰ = FT(1), αᵀ = FT(0), κʰ = FT(0), @@ -141,13 +147,13 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ()) model_2D = ShallowWaterModel( param_set, prob_2D, + coupling, ShallowWater.ConstantViscosity{FT}(model_3D.νʰ), nothing, FT(1), ) dt_fast = 300 - dt_slow = 300 nout = ceil(Int64, tout / dt_slow) dt_slow = tout / nout @@ -180,7 +186,7 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ()) dg_2D = DGModel( model_2D, grid_2D, - RusanovNumericalFlux(), + CentralNumericalFluxFirstOrder(), CentralNumericalFluxSecondOrder(), CentralNumericalFluxGradient(), ) @@ -190,8 +196,7 @@ function run_hydrostatic_spindown(; coupled = true, refDat = ()) lsrk_3D = LSRK54CarpenterKennedy(dg_3D, Q_3D, dt = dt_slow, t0 = 0) lsrk_2D = LSRK54CarpenterKennedy(dg_2D, Q_2D, dt = dt_fast, t0 = 0) - - odesolver = SplitExplicitSolver(lsrk_3D, lsrk_2D; coupled = coupled) + odesolver = SplitExplicitSolver(lsrk_3D, lsrk_2D) step = [0, 0] cbvector = make_callbacks( @@ -312,6 +317,7 @@ function make_callbacks( (Q_slow, "3D state"), (dg_slow.state_auxiliary, "3D aux"), (Q_fast, "2D state"), + (dg_fast.state_auxiliary, "2D aux"), ], nout; prec = 12, @@ -327,7 +333,7 @@ FT = Float64 vtkpath = "vtk_split" const timeend = 24 * 3600 # s -const tout = 2 * 3600 # s +const tout = 3 * 3600 # s # const timeend = 1200 # s # const tout = 300 # s @@ -350,5 +356,52 @@ const cᶻ = 0 @testset "$(@__FILE__)" begin include("../refvals/hydrostatic_spindown_refvals.jl") - run_hydrostatic_spindown(coupled = false, refDat = refVals.uncoupled) + @testset "Multi-rate" begin + + @testset "Fully Coupled" begin + @testset "Δt = 30 mins" begin + run_hydrostatic_spindown( + coupling = Coupled(), + dt_slow = 30 * 60, + refDat = refVals.thirty_minutes, + ) + end + + @testset "Δt = 60 mins" begin + run_hydrostatic_spindown( + coupling = Coupled(), + dt_slow = 60 * 60, + refDat = refVals.sixty_minutes, + ) + end + + @testset "Δt = 90 mins" begin + run_hydrostatic_spindown( + coupling = Coupled(), + dt_slow = 90 * 60, + refDat = refVals.ninety_minutes, + ) + end + end + end + + if ClimateMachine.Settings.integration_testing + @testset "Single-Rate" begin + @testset "Not Coupled" begin + run_hydrostatic_spindown( + coupling = Uncoupled(), + dt_slow = 300, + refDat = refVals.not_coupled, + ) + end + + @testset "Fully Coupled" begin + run_hydrostatic_spindown( + coupling = Coupled(), + dt_slow = 300, + refDat = refVals.fully_coupled, + ) + end + end + end end diff --git a/test/Ocean/SplitExplicit/test_vertical_integral_model.jl b/test/Ocean/SplitExplicit/test_vertical_integral_model.jl index 9d51ea69dfd..483e8c9d80c 100644 --- a/test/Ocean/SplitExplicit/test_vertical_integral_model.jl +++ b/test/Ocean/SplitExplicit/test_vertical_integral_model.jl @@ -7,6 +7,7 @@ using ClimateMachine.ODESolvers using ClimateMachine.Mesh.Filters using ClimateMachine.VariableTemplates using ClimateMachine.Mesh.Grids: polynomialorder +using ClimateMachine.Ocean using ClimateMachine.Ocean.HydrostaticBoussinesq using ClimateMachine.Ocean.ShallowWater using ClimateMachine.Ocean.SplitExplicit: VerticalIntegralModel @@ -131,6 +132,7 @@ function test_vertical_integral_model(time; refDat = ()) model_2D = ShallowWaterModel( param_set, prob_2D, + Uncoupled(), ShallowWater.ConstantViscosity{FT}(model_3D.νʰ), nothing, FT(1),