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))