From 6f34be46bd99c27d5694d5a252425eac1b10fe8e Mon Sep 17 00:00:00 2001 From: Brandon Allen Date: Thu, 25 Jun 2020 13:59:47 -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 --- .../ODESolvers/SplitExplicitMethod.jl | 125 ++++------ .../HydrostaticBoussinesqModel.jl | 99 +++++--- src/Ocean/Ocean.jl | 8 + src/Ocean/ShallowWater/ShallowWaterModel.jl | 38 ++- src/Ocean/SplitExplicit/Communication.jl | 46 ++-- .../HydrostaticBoussinesqCoupling.jl | 81 +++++++ .../SplitExplicit/ShallowWaterCoupling.jl | 35 +++ src/Ocean/SplitExplicit/SplitExplicitModel.jl | 42 ++++ .../ShallowWater/2D_hydrostatic_spindown.jl | 2 + test/Ocean/ShallowWater/GyreDriver.jl | 9 +- test/Ocean/ShallowWater/GyreInABox.jl | 1 + .../SplitExplicit/hydrostatic_spindown.jl | 71 +++--- .../test_vertical_integral_model.jl | 2 +- .../refvals/hydrostatic_spindown_refvals.jl | 220 +++++++++++++++++- 14 files changed, 593 insertions(+), 186 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 0520c8f7f2e..78b253d9129 100644 --- a/src/Numerics/ODESolvers/SplitExplicitMethod.jl +++ b/src/Numerics/ODESolvers/SplitExplicitMethod.jl @@ -27,7 +27,7 @@ of the model at a faster rate than the full model. This results in a first- order time stepper. """ SplitExplicitSolver -mutable struct SplitExplicitSolver{SS, FS, RT, BT, MSA} <: AbstractODESolver +mutable struct SplitExplicitSolver{SS, FS, RT, MSA} <: AbstractODESolver "slow solver" slow_solver::SS "fast solver" @@ -36,10 +36,6 @@ mutable struct SplitExplicitSolver{SS, FS, RT, BT, MSA} <: AbstractODESolver dt::RT "time" t::RT - "flag to drive fast model by slow model" - coupled::BT - "flag to reconcile slow model with fast model" - reconcile::BT "storage for transfer tendency" dQ2fast::MSA @@ -49,25 +45,19 @@ mutable struct SplitExplicitSolver{SS, FS, RT, BT, MSA} <: AbstractODESolver Q = nothing; dt = getdt(slow_solver), t0 = slow_solver.t, - coupled = true, - reconcile = true, ) where {AT <: AbstractArray} SS = typeof(slow_solver) FS = typeof(fast_solver) RT = real(eltype(slow_solver.dQ)) - BT = typeof(coupled) - dQ2fast = similar(slow_solver.dQ) dQ2fast .= -0 MSA = typeof(dQ2fast) - return new{SS, FS, RT, BT, MSA}( + return new{SS, FS, RT, MSA}( slow_solver, fast_solver, RT(dt), RT(t0), - coupled, - reconcile, dQ2fast, ) end @@ -102,54 +92,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, - split.reconcile, - ) - 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] @@ -167,31 +131,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, - split.reconcile, + 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..48fdda963d2 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! @@ -18,6 +19,8 @@ using ...DGMethods.NumericalFluxes: RusanovNumericalFlux using ...BalanceLaws using ...BalanceLaws: number_states +import ..Ocean: coriolis_force! +import ...DGMethods.NumericalFluxes: update_penalty! import ...BalanceLaws: vars_state, init_state_prognostic!, @@ -38,6 +41,7 @@ import ...BalanceLaws: reverse_integral_load_auxiliary_state!, reverse_integral_set_auxiliary_state! + ×(a::SVector, b::SVector) = StaticArrays.cross(a, b) ⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b) ⊗(a::SVector, b::SVector) = a * b' @@ -67,9 +71,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 @@ -84,6 +89,7 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw function HydrostaticBoussinesqModel{FT}( param_set::PS, problem; + coupling = Uncoupled(), ρₒ = FT(1000), # kg / m^3 cʰ = FT(0), # m/s cᶻ = FT(0), # m/s @@ -96,9 +102,10 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw fₒ = FT(1e-4), # Hz β = FT(1e-11), # Hz / m ) where {FT <: AbstractFloat, PS} - return new{PS, typeof(problem), FT}( + return new{typeof(coupling), PS, typeof(problem), FT}( param_set, problem, + coupling, ρₒ, cʰ, cᶻ, @@ -189,6 +196,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 +215,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 +274,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 +282,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 +450,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 +546,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 +572,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 +632,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) diff --git a/src/Ocean/Ocean.jl b/src/Ocean/Ocean.jl index 9dba464b5f5..412063993e1 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_force! 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..569d3e2332b 100644 --- a/src/Ocean/ShallowWater/ShallowWaterModel.jl +++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl @@ -6,12 +6,14 @@ using StaticArrays using LinearAlgebra: dot, Diagonal using CLIMAParameters.Planet: grav +using ..Ocean using ...VariableTemplates using ...Mesh.Geometry using ...DGMethods using ...DGMethods.NumericalFluxes using ...BalanceLaws +import ..Ocean: coriolis_force! import ...DGMethods.NumericalFluxes: update_penalty! import ...BalanceLaws: vars_state, @@ -25,6 +27,7 @@ import ...BalanceLaws: wavespeed, boundary_state! + ×(a::SVector, b::SVector) = StaticArrays.cross(a, b) ⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b) ⊗(a::SVector, b::SVector) = a * b' @@ -42,9 +45,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 +160,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,13 +191,15 @@ function flux_second_order!( α::Vars, t::Real, ) - flux_second_order!(m.turbulence, G, q, σ, α, t) + flux_second_order!(m, m.turbulence, m.coupling, G, q, σ, α, t) end -flux_second_order!(::LinearDrag, _...) = nothing +flux_second_order!(::SWModel, ::LinearDrag, _...) = nothing @inline function flux_second_order!( + ::SWModel, ::ConstantViscosity, + ::Uncoupled, G::Grad, q::Vars, σ::Vars, @@ -210,18 +217,29 @@ end m::SWModel{P}, S::Vars, q::Vars, - diffusive::Vars, + d::Vars, α::Vars, t::Real, direction, ) where {P} - S.U += α.τ + coriolis_force!(m, m.coupling, S, q, α, t) + forcing_term!(m, m.coupling, S, q, α, t) + linear_drag!(m.turbulence, S, q, α, t) + + return nothing +end - f = α.f - U, V = q.U +@inline function coriolis_force!(::SWModel, ::Uncoupled, S, Q, A, t) + # f × u + f = A.f + U, V = Q.U S.U -= @SVector [-f * V, f * U] - 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 diff --git a/src/Ocean/SplitExplicit/Communication.jl b/src/Ocean/SplitExplicit/Communication.jl index 1e088c8973a..49ff2305fb4 100644 --- a/src/Ocean/SplitExplicit/Communication.jl +++ b/src/Ocean/SplitExplicit/Communication.jl @@ -1,32 +1,25 @@ -import ...BalanceLaws: - initialize_states!, - tendency_from_slow_to_fast!, - cummulate_fast_solution!, - reconcile_from_fast_to_slow! - @inline function initialize_states!( - baroclinic::HydrostaticBoussinesqModel, - barotropic::ShallowWaterModel, + baroclinic::HBModel{C}, + barotropic::SWModel{C}, model_bc, model_bt, state_bc, state_bt, -) +) where {C <: Coupled} model_bc.state_auxiliary.ΔGᵘ .= -0 return nothing end @inline function tendency_from_slow_to_fast!( - baroclinic::HydrostaticBoussinesqModel, - barotropic::ShallowWaterModel, + baroclinic::HBModel{C}, + barotropic::SWModel{C}, model_bc, model_bt, state_bc, state_bt, forcing_tendency, - reconcile, -) +) where {C <: Coupled} FT = eltype(state_bc) Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(model_bc) @@ -60,35 +53,37 @@ end Gᵁ = @view data_bt[:, index_Gᵁ, 1:nrealelemh] Gᵁ .= ∫du + ### get top value (=integral over full depth) of ∫du + ∫du = @view data_int[:, end:end, index_∫du, end:end, 1:nrealelemh] + ### save vertically averaged tendency to remove from 3D tendency ### need to reshape for the broadcast ΔGᵘ = @view data_bc[:, :, index_ΔGᵘ, :, 1:nrealelemh] - ΔGᵘ .-= reshape(∫du, Nq^2, 1, 2, 1, nrealelemh) / baroclinic.problem.H + ΔGᵘ .-= ∫du / baroclinic.problem.H return nothing end @inline function cummulate_fast_solution!( - baroclinic::HydrostaticBoussinesqModel, - barotropic::ShallowWaterModel, + baroclinic::HBModel{C}, + barotropic::SWModel{C}, model_bt, state_bt, fast_time, fast_dt, substep, -) +) where {C <: Coupled} return nothing end @inline function reconcile_from_fast_to_slow!( - baroclinic::HydrostaticBoussinesqModel, - barotropic::ShallowWaterModel, + baroclinic::HBModel{C}, + barotropic::SWModel{C}, model_bc, model_bt, state_bc, state_bt, - reconcile, -) +) where {C <: Coupled} FT = eltype(state_bc) Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(model_bc) @@ -133,14 +128,17 @@ end ### copy the 2D contribution down the 3D solution ### need to reshape for the broadcast + data_bt_aux = reshape(data_bt_aux, Nq^2, 1, num_aux_bt, 1, nelemh) + Δu = @view data_bt_aux[:, :, index_Δu, :, 1:nrealelemh] u = @view data_bc_state[:, :, index_u, :, 1:nrealelemh] - u .+= reshape(Δu, Nq^2, 1, 2, 1, nrealelemh) + u .+= Δu ### copy η from barotropic mode to baroclinic mode ### need to reshape for the broadcast + data_bt_state = reshape(data_bt_state, Nq^2, 1, num_state_bt, 1, nelemh) + η_2D = @view data_bt_state[:, :, index_η_2D, :, 1:nrealelemh] η_3D = @view data_bc_state[:, :, index_η_3D, :, 1:nrealelemh] - η_2D = @view data_bt_state[:, index_η_2D, 1:nrealelemh] - η_3D .= reshape(η_2D, Nq^2, 1, 1, 1, nrealelemh) + η_3D .= η_2D return nothing end 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..314628fd4fe --- /dev/null +++ b/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl @@ -0,0 +1,35 @@ +using ..ShallowWater +using ..ShallowWater: ConstantViscosity + +import ...BalanceLaws: flux_second_order! +import ..ShallowWater: forcing_term! + +@inline function flux_second_order!( + ::SWModel, + ::ConstantViscosity, + ::Coupled, + G::Grad, + q::Vars, + σ::Vars, + α::Vars, + t::Real, +) + G.U += σ.ν∇U + + return nothing +end + +@inline function coriolis_force!(::SWModel, ::Coupled, S, Q, A, t) + # f × u + f = A.f + U, V = Q.U + S.U -= @SVector [-f * V, f * U] + + return nothing +end + +@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 4bb3e67fe27..9efeac21464 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 @@ -11,7 +12,48 @@ using ...Mesh.Geometry using ...DGMethods using ...BalanceLaws +import ..Ocean: coriolis_force! +import ...BalanceLaws: + initialize_states!, + tendency_from_slow_to_fast!, + 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 94ecd12334f..0164a509598 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 @@ -86,8 +87,8 @@ function shallow_init_aux!(m::ShallowWaterModel, p::GyreInABox, A, geom) end function run_hydrostatic_spindown(; - coupled = true, - reconcile = true, + coupling = Uncoupled(), + dt_slow = 300, refDat = (), ) mpicomm = MPI.COMM_WORLD @@ -134,6 +135,7 @@ function run_hydrostatic_spindown(; model_3D = HydrostaticBoussinesqModel{FT}( param_set, prob_3D; + coupling = coupling, cʰ = FT(1), αᵀ = FT(0), κʰ = FT(0), @@ -145,13 +147,13 @@ function run_hydrostatic_spindown(; 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 @@ -194,8 +196,7 @@ function run_hydrostatic_spindown(; 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( @@ -314,8 +315,9 @@ function make_callbacks( cbcs_dg = ClimateMachine.StateCheck.sccreate( [ (Q_slow, "3D state"), - # dg_slow.state_auxiliary, "3D aux"), + # (dg_slow.state_auxiliary, "3D aux"), (Q_fast, "2D state"), + # (dg_fast.state_auxiliary, "2D aux"), ], nout; prec = 12, @@ -331,9 +333,9 @@ FT = Float64 vtkpath = "vtk_split" const timeend = 24 * 3600 # s -const tout = 2 * 3600 # s -const timeend = 1200 # s -const tout = 300 # s +const tout = 3 * 3600 # s +# const timeend = 1200 # s +# const tout = 300 # s const N = 4 const Nˣ = 5 @@ -354,43 +356,21 @@ const cᶻ = 0 @testset "$(@__FILE__)" begin include("../refvals/hydrostatic_spindown_refvals.jl") - #= - @testset "Single-Rate" begin - @testset "Not Coupled" begin - run_hydrostatic_spindown( - coupling = NotCoupled(), - dt_slow = 300, - refDat = refVals.not_coupled, - ) - end - - # run_hydrostatic_spindown(coupling = PartiallyCoupled()) - - @testset "Fully Coupled" begin - run_hydrostatic_spindown( - coupling = FullyCoupled(), - dt_slow = 300, - refDat = refVals.fully_coupled, - ) - end - end - =# - @testset "Multi-rate" begin # run_hydrostatic_spindown(coupling = PartiallyCoupled()) @testset "Fully Coupled" begin @testset "Δt = 30 mins" begin run_hydrostatic_spindown( - coupling = FullyCoupled(), + coupling = Coupled(), dt_slow = 30 * 60, refDat = refVals.thirty_minutes, ) end - #= + @testset "Δt = 60 mins" begin run_hydrostatic_spindown( - coupling = FullyCoupled(), + coupling = Coupled(), dt_slow = 60 * 60, refDat = refVals.sixty_minutes, ) @@ -398,12 +378,31 @@ const cᶻ = 0 @testset "Δt = 90 mins" begin run_hydrostatic_spindown( - coupling = FullyCoupled(), + 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 ba71500e72a..d5c07492faf 100644 --- a/test/Ocean/SplitExplicit/test_vertical_integral_model.jl +++ b/test/Ocean/SplitExplicit/test_vertical_integral_model.jl @@ -132,7 +132,7 @@ function test_vertical_integral_model(time; refDat = ()) model_2D = ShallowWaterModel( param_set, prob_2D, - NotCoupled(), + Uncoupled(), ShallowWater.ConstantViscosity{FT}(model_3D.νʰ), nothing, FT(1), diff --git a/test/Ocean/refvals/hydrostatic_spindown_refvals.jl b/test/Ocean/refvals/hydrostatic_spindown_refvals.jl index 00f32a3bc72..20a39c2cf90 100644 --- a/test/Ocean/refvals/hydrostatic_spindown_refvals.jl +++ b/test/Ocean/refvals/hydrostatic_spindown_refvals.jl @@ -3,20 +3,35 @@ # [ : : : : : : ], # ] parr = [ - ["Q", "u[1]", 12, 12, 0, 12], - ["Q", "u[2]", 0, 0, 0, 0], - ["Q", :η, 12, 12, 0, 12], - ["Q", :θ, 15, 15, 15, 15], - # ["s_aux", :y, 15, 15, 15, 15], - # ["s_aux", :w, 12, 12, 0, 12], - # ["s_aux", :pkin, 15, 15, 15, 15], - # ["s_aux", :wz0, 12, 12, 0, 12], + ["3D state", "u[1]", 12, 12, 0, 12], + ["3D state", "u[2]", 0, 0, 0, 0], + ["3D state", :η, 12, 12, 0, 12], + ["3D state", :θ, 15, 15, 15, 15], + #= + ["3D aux", :y, 15, 15, 15, 15], + ["3D aux", :w, 12, 12, 0, 12], + ["3D aux", :pkin, 15, 15, 15, 15], + ["3D aux", :wz0, 12, 12, 0, 12], + ["3D aux", "uᵈ[1]", 12, 12, 0, 12], + ["3D aux", "uᵈ[2]", 12, 12, 0, 12], + ["3D aux", "ΔGᵘ[1]", 12, 12, 0, 12], + ["3D aux", "ΔGᵘ[2]", 12, 12, 0, 12], + =# ["2D state", :η, 12, 12, 0, 12], ["2D state", "U[1]", 12, 12, 0, 12], ["2D state", "U[2]", 0, 0, 0, 0], + #= + ["2D aux", :f, 15, 15, 15, 15], + ["2D aux", "τ[1]", 15, 15, 15, 15], + ["2D aux", "τ[2]", 15, 15, 15, 15], + ["2D aux", "Gᵁ[1]", 12, 12, 0, 12], + ["2D aux", "Gᵁ[2]", 12, 12, 0, 12], + ["2D aux", "Δu[1]", 12, 12, 0, 12], + ["2D aux", "Δu[2]", 12, 12, 0, 12], + =# ] -uncoupled = [ +not_coupled = [ [ "3D state", "u[1]", @@ -109,4 +124,189 @@ uncoupled = [ ], ] -refVals = (uncoupled = (uncoupled, parr),) +fully_coupled = [ + [ + "3D state", + "u[1]", + -9.58660726483847924762e-01, + 9.58615148375367986944e-01, + -3.89899325474971198514e-06, + 4.45401015112516618366e-01, + ], + [ + "3D state", + "u[2]", + -5.34954449329680479380e-05, + 9.56394271461731787596e-05, + 1.41658660512563957739e-07, + 1.03060542771323416372e-05, + ], + [ + "3D state", + :η, + -8.53000730930055017787e-01, + 8.52896671352524671228e-01, + 1.16182894068970207099e-06, + 6.02991581555750166821e-01, + ], + [ + "3D state", + :θ, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + ], + #= + [ + "3D aux", + :y, + 0.00000000000000000000e+00, + 1.00000000000000011642e+06, + 5.00000000000000000000e+05, + 2.92775877460665535182e+05, + ], + [ + "3D aux", + :w, + -3.53444692332414609707e-04, + 3.53393768653404403979e-04, + -2.62684474750210507396e-09, + 1.76703229531742088583e-04, + ], + [ + "3D aux", + :pkin, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + ], + [ + "3D aux", + :wz0, + -2.75669561310449110355e-16, + 2.61116688526775137768e-16, + 1.59449559147079600665e-18, + 3.68498533954047109486e-17, + ], + [ + "3D aux", + "uᵈ[1]", + -8.79815521725190063940e-01, + 8.79860306849815310137e-01, + -3.36235206981427825909e-08, + 4.41872889723482875635e-01, + ], + [ + "3D aux", + "uᵈ[2]", + -7.64485889244402041197e-08, + 1.07475092686394092540e-07, + 1.97594920215074170770e-12, + 1.00383435739831354977e-08, + ], + [ + "3D aux", + "ΔGᵘ[1]", + -8.90270808298916137079e-09, + 8.90270808302469374118e-09, + -3.03471576210481703598e-21, + 1.85087492582869203479e-09, + ], + [ + "3D aux", + "ΔGᵘ[2]", + -3.15722602438961622861e-09, + 5.35717978250396107350e-09, + -2.15993401549846601979e-25, + 7.02233250473352312290e-10, + ], + =# + [ + "2D state", + :η, + -8.53000730930055017787e-01, + 8.52896671352524671228e-01, + 1.16182894073002528620e-06, + 6.03462484854037639614e-01, + ], + [ + "2D state", + "U[1]", + -3.15731740631152391074e+01, + 3.15623403633881167707e+01, + -1.54618387206766983236e-03, + 2.24243954053573801843e+01, + ], + [ + "2D state", + "U[2]", + -2.13769342942716975009e-02, + 3.82134429332010097657e-02, + 5.66626758280705091560e-05, + 4.12563894411444639226e-03, + ], + #= + [ + "2D aux", + :f, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + ], + [ + "2D aux", + "τ[1]", + -0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + ], + [ + "2D aux", + "τ[2]", + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + 0.00000000000000000000e+00, + ], + [ + "2D aux", + "Gᵁ[1]", + -3.56108323320987754941e-06, + 3.56108323319566433656e-06, + 1.21388901534735797970e-18, + 7.40928142985258653790e-07, + ], + [ + "2D aux", + "Gᵁ[2]", + -2.14287191300158429705e-06, + 1.26289040975584651791e-06, + 9.16489648929153005052e-23, + 2.81112662425186374727e-07, + ], + [ + "2D aux", + "Δu[1]", + -6.54477648201510250718e-04, + 6.54167531208154056505e-04, + 1.02401396398448558772e-07, + 4.64940295021459143454e-04, + ], + [ + "2D aux", + "Δu[2]", + -2.30864769013247066032e-06, + 4.06331319646622782525e-06, + 9.84808153015545305074e-10, + 4.00443091670026613147e-07, + ], + =# +] + + +refVals = + (not_coupled = (not_coupled, parr), fully_coupled = (fully_coupled, parr))