diff --git a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl index a4e4da8360f..732279f5a94 100644 --- a/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl +++ b/src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl @@ -14,6 +14,7 @@ using ...Mesh.Geometry using ...DGMethods: DGModel, copy_stack_field_down! using ...DGMethods.NumericalFluxes using ...BalanceLaws +using ..Ocean import ...DGMethods.NumericalFluxes: update_penalty! import ...BalanceLaws: @@ -40,6 +41,7 @@ import ...BalanceLaws: reverse_indefinite_stack_integral!, reverse_integral_load_auxiliary_state!, reverse_integral_set_auxiliary_state! +import ..Ocean: coriolis_force! ×(a::SVector, b::SVector) = StaticArrays.cross(a, b) ⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b) @@ -70,9 +72,10 @@ fₒ = first coriolis parameter (constant term) HydrostaticBoussinesqModel(problem) """ -struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw +struct HydrostaticBoussinesqModel{PS, P, C, T} <: BalanceLaw param_set::PS problem::P + coupling::C ρₒ::T cʰ::T cᶻ::T @@ -86,6 +89,7 @@ struct HydrostaticBoussinesqModel{PS, P, T} <: BalanceLaw function HydrostaticBoussinesqModel{FT}( param_set::PS, problem; + coupling = NotCoupled(), ρₒ = FT(1000), # kg / m^3 cʰ = FT(0), # m/s cᶻ = FT(0), # m/s @@ -97,9 +101,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{PS, typeof(problem), typeof(coupling), FT}( param_set, problem, + coupling, ρₒ, cʰ, cᶻ, @@ -189,6 +194,7 @@ these are just copies in our model function vars_state_gradient(m::HBModel, T) @vars begin ∇u::SVector{2, T} + ∇uᵈ::SVector{2, T} ∇θ::T end end @@ -207,9 +213,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, + ::NotCoupled, + G, + Q, + A, + t, +) + G.∇u = Q.u + return nothing end @@ -248,8 +268,7 @@ this computation is done pointwise at each nodal point A::Vars, t, ) - ν = 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.∇θ @@ -257,6 +276,21 @@ this computation is done pointwise at each nodal point return nothing end +@inline function velocity_gradient_flux!( + m::HBModel, + ::NotCoupled, + D, + G, + Q, + A, + t, +) + ν = viscosity_tensor(m) + D.ν∇u = -ν * G.∇u + + return nothing +end + """ viscosity_tensor(::HBModel) @@ -419,38 +453,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, ::NotCoupled, 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) @@ -509,22 +549,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, ::NotCoupled, 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 @@ -532,7 +575,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) @@ -592,9 +635,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, ::NotCoupled, _...) = nothing + """ update_auxiliary_state_gradient!(::HBModel) diff --git a/src/Ocean/Ocean.jl b/src/Ocean/Ocean.jl index 9dba464b5f5..62d3c176031 100644 --- a/src/Ocean/Ocean.jl +++ b/src/Ocean/Ocean.jl @@ -1,5 +1,14 @@ module Ocean +export AbstractOceanCoupling, NotCoupled, PartiallyCoupled, FullyCoupled + +abstract type AbstractOceanCoupling end +struct NotCoupled <: AbstractOceanCoupling end +struct PartiallyCoupled <: AbstractOceanCoupling end +struct FullyCoupled <: 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 fc2f0ecab8c..63be50facec 100644 --- a/src/Ocean/ShallowWater/ShallowWaterModel.jl +++ b/src/Ocean/ShallowWater/ShallowWaterModel.jl @@ -11,6 +11,7 @@ using ...Mesh.Geometry using ...DGMethods using ...DGMethods.NumericalFluxes using ...BalanceLaws +using ..Ocean import ...DGMethods.NumericalFluxes: update_penalty! import ...BalanceLaws: @@ -27,6 +28,7 @@ import ...BalanceLaws: source!, wavespeed, boundary_state! +import ..Ocean: coriolis_force! ×(a::SVector, b::SVector) = StaticArrays.cross(a, b) ⋅(a::SVector, b::SVector) = StaticArrays.dot(a, b) @@ -45,9 +47,10 @@ end abstract type AdvectionTerm end struct NonLinearAdvection <: AdvectionTerm end -struct ShallowWaterModel{PS, P, T, A, S} <: BalanceLaw +struct ShallowWaterModel{PS, P, T, A, C, S} <: BalanceLaw param_set::PS problem::P + coupling::C turbulence::T advection::A c::S @@ -159,12 +162,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, @@ -189,13 +193,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, + ::NotCoupled, G::Grad, q::Vars, σ::Vars, @@ -213,18 +219,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, ::NotCoupled, 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, ::NotCoupled, S, Q, A, t) + S.U += A.τ return nothing end diff --git a/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl b/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl new file mode 100644 index 00000000000..4ebbd86cf0a --- /dev/null +++ b/src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl @@ -0,0 +1,141 @@ +using ..HydrostaticBoussinesq +using ...DGMethods + +import ...BalanceLaws +import ..HydrostaticBoussinesq: + viscosity_tensor, + coriolis_parameter, + velocity_gradient_argument!, + velocity_gradient_flux!, + hydrostatic_pressure!, + compute_flow_deviation! + +HBModel = HydrostaticBoussinesqModel + +@inline function velocity_gradient_argument!( + m::HBModel, + ::PartiallyCoupled, + G, + Q, + A, + t, +) + G.∇u = Q.u + + return nothing +end + +@inline function velocity_gradient_argument!( + m::HBModel, + ::FullyCoupled, + G, + Q, + A, + t, +) + G.∇u = Q.u + G.∇uᵈ = A.uᵈ + + return nothing +end + +@inline function velocity_gradient_flux!( + m::HBModel, + ::PartiallyCoupled, + D, + G, + Q, + A, + t, +) + ν = viscosity_tensor(m) + D.ν∇u = -ν * G.∇u + + return nothing +end + +@inline function velocity_gradient_flux!( + m::HBModel, + ::FullyCoupled, + 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 function hydrostatic_pressure!( + m::HBModel, + ::PartiallyCoupled, + F, + Q, + A, + t, +) + η = Q.η + Iʰ = @SMatrix [ + 1 -0 + -0 1 + -0 -0 + ] + + F.u += grav(m.param_set) * η * Iʰ + + return nothing +end + +@inline hydrostatic_pressure!(::HBModel, ::FullyCoupled, _...) = nothing + +@inline function coriolis_force!(m::HBModel, ::PartiallyCoupled, 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 + +@inline function coriolis_force!(m::HBModel, ::FullyCoupled, 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 + +@inline compute_flow_deviation!(dg, ::HBModel, ::PartiallyCoupled, _...) = + nothing + +# Compute Horizontal Flow deviation from vertical mean +@inline function compute_flow_deviation!(dg, m::HBModel, ::FullyCoupled, Q, t) + Nq, Nqk, _, _, nelemv, _, nelemh, _ = basic_grid_info(dg) + + flow_model = dg.modeldata.flow_model + update_auxiliary_state!(flow_model, flow_model.balance_law, Q, 0) + + ## get top value (=integral over full depth) + ∫u = flow_model.state_auxiliary.∫u + boxy_∫u = reshape(∫u, Nq^2, Nqk, 2, nelemv, nelemh) + flat_∫u = @view boxy_∫u[:, end, :, end, :] + boxy_uᵇ = reshape(flat_∫u, Nq^2, 1, 2, 1, nelemh) + + ## make a copy of horizontal velocity + ## and remove vertical mean velocity + uᵈ = dg.state_auxiliary.uᵈ + uᵈ .= Q.u + boxy_uᵈ = reshape(uᵈ, Nq^2, Nqk, 2, nelemv, nelemh) + boxy_uᵈ .-= boxy_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..0b9e1b17480 --- /dev/null +++ b/src/Ocean/SplitExplicit/ShallowWaterCoupling.jl @@ -0,0 +1,58 @@ +using ..ShallowWater +using ..ShallowWater: ConstantViscosity + +import ...BalanceLaws: flux_second_order! +import ..ShallowWater: forcing_term! + +flux_second_order!( + ::ShallowWaterModel, + ::ConstantViscosity, + ::PartiallyCoupled, + _..., +) = nothing + +@inline function flux_second_order!( + ::ShallowWaterModel, + ::ConstantViscosity, + ::FullyCoupled, + G::Grad, + q::Vars, + σ::Vars, + α::Vars, + t::Real, +) + G.U += σ.ν∇U + + return nothing +end + +coriolis_force!(::ShallowWaterModel, ::PartiallyCoupled, _...) = nothing + +@inline function coriolis_force!( + ::ShallowWaterModel, + ::FullyCoupled, + 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!( + ::ShallowWaterModel, + ::Union{PartiallyCoupled, FullyCoupled}, + 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 bac8641ab31..606a5c35f1e 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,6 +12,7 @@ using ...Mesh.Geometry using ...DGMethods using ...BalanceLaws +import ..Ocean: coriolis_force! import ...BalanceLaws: vars_state_conservative, vars_state_gradient, @@ -26,6 +28,7 @@ import ...BalanceLaws: include("TendencyIntegralModel.jl") include("FlowIntegralModel.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 4bd776b23fe..3fb19db7925 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, + NotCoupled(), ShallowWater.ConstantViscosity{FT}(5e3), nothing, FT(1), diff --git a/test/Ocean/ShallowWater/GyreDriver.jl b/test/Ocean/ShallowWater/GyreDriver.jl index 96b20d22560..5b7e799d2bc 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, + NotCoupled(), + turbulence, + advection, + c, + ) end function shallow_init_state!( diff --git a/test/Ocean/ShallowWater/GyreInABox.jl b/test/Ocean/ShallowWater/GyreInABox.jl index fc3716922dc..8b5cdbefce0 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 ede2172d799..ae33156af74 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: @@ -87,8 +88,8 @@ function shallow_init_aux!(m::ShallowWaterModel, p::GyreInABox, A, geom) end function run_hydrostatic_spindown(; - coupled = true, - reconcile = true, + coupling = NotCoupled(), + dt_slow = 300, refDat = (), ) mpicomm = MPI.COMM_WORLD @@ -135,6 +136,7 @@ function run_hydrostatic_spindown(; model_3D = HydrostaticBoussinesqModel{FT}( param_set, prob_3D; + coupling = coupling, cʰ = FT(1), αᵀ = FT(0), κʰ = FT(0), @@ -146,13 +148,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 @@ -205,7 +207,23 @@ 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) + if typeof(coupling) == NotCoupled + coupled = false + reconcile = false + elseif typeof(coupling) == PartiallyCoupled + coupled = true + reconcile = false + elseif typeof(coupling) == FullyCoupled + coupled = true + reconcile = true + end + + odesolver = SplitExplicitSolver( + lsrk_3D, + lsrk_2D; + coupled = coupled, + reconcile = reconcile, + ) step = [0, 0] cbvector = make_callbacks( @@ -324,8 +342,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, @@ -341,9 +360,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 @@ -364,12 +383,51 @@ const cᶻ = 0 @testset "$(@__FILE__)" begin include("../refvals/hydrostatic_spindown_refvals.jl") - #= - run_hydrostatic_spindown( - coupled = false, - reconcile = false, - refDat = refVals.uncoupled, - ) - =# - run_hydrostatic_spindown(coupled = true, reconcile = false) + @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(), + dt_slow = 30 * 60, + # refDat = refVals.fully_coupled, + ) + end + @testset "Δt = 60 mins" begin + run_hydrostatic_spindown( + coupling = FullyCoupled(), + dt_slow = 60 * 60, + # refDat = refVals.fully_coupled, + ) + end + @testset "Δt = 90 mins" begin + run_hydrostatic_spindown( + coupling = FullyCoupled(), + dt_slow = 90 * 60, + # refDat = refVals.fully_coupled, + ) + end + end + end end 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))