Merge #1515

1515: Add rotating spindown test to handling of Coriolis force r=blallen a=blallen

# Description

Another spinoff from #1482. Depends on #1509. Adds a `rotation` field to the `SimpleBox` problem to turn on and off the Coriolis force in the hydrostatic spindown tests. Need to add state check values after the pretty-tables output is fixed.



Co-authored-by: blallen <[email protected]>

Author: bors[bot]

.buildkite/pipeline.yml

@@ -381,6 +381,15 @@ steps:
       queue: central
       slurm_ntasks: 1
 
+  - label: "cpu_test_coriolis"
+    key: "cpu_test_coriolis"
+    command:
+      - "mpiexec julia --color=yes --project test/Ocean/SplitExplicit/test_coriolis.jl "
+    agents:
+      config: cpu
+      queue: central
+      slurm_ntasks: 1
+
   - label: "cpu_KM_saturation_adjustment"
     key: "cpu_KM_saturation_adjustment"
     command:
@@ -754,6 +763,16 @@ steps:
       slurm_ntasks: 1
       slurm_gres: "gpu:1"
 
+  - label: "gpu_test_coriolis"
+    key: "gpu_test_coriolis"
+    command:
+      - "mpiexec julia --color=yes --project test/Ocean/SplitExplicit/test_coriolis.jl "
+    agents:
+      config: gpu
+      queue: central
+      slurm_ntasks: 1
+      slurm_gres: "gpu:1"
+
   - label: "gpu_KM_saturation_adjustment"
     key: "gpu_KM_saturation_adjustment"
     command:

src/Ocean/HydrostaticBoussinesq/HydrostaticBoussinesqModel.jl

@@ -607,24 +607,13 @@ end
 
 @inline function coriolis_force!(m::HBModel, ::Uncoupled, S, Q, A, t)
     # f × u
-    f = coriolis_parameter(m, A.y)
+    f = coriolis_parameter(m, m.problem, A.y)
     u, v = Q.u # Horizontal components of velocity
     S.u -= @SVector [-f * v, f * u]
 
     return nothing
 end
 
-"""
-    coriolis_parameter(::HBModel)
-
-northern hemisphere coriolis
-
-# Arguments
-- `m`: model object to dispatch on and get coriolis parameters
-- `y`: y-coordinate in the box
-"""
-@inline coriolis_parameter(m::HBModel, y) = m.fₒ + m.β * y
-
 """
     wavespeed(::HBModel)
 

src/Ocean/OceanProblems/SimpleBoxProblem.jl

@@ -1,6 +1,6 @@
 module OceanProblems
 
-export SimpleBox, HomogeneousBox, OceanGyre
+export SimpleBox, Fixed, Rotating, HomogeneousBox, OceanGyre
 
 using StaticArrays
 using CLIMAParameters.Planet: grav
@@ -23,10 +23,6 @@ SWModel = ShallowWaterModel
 
 abstract type AbstractOceanProblem <: AbstractProblem end
 
-############################
-# Basic box problem        #
-# Set up dimensions of box #
-############################
 abstract type AbstractSimpleBoxProblem <: AbstractOceanProblem end
 
 """
@@ -53,8 +49,6 @@ function ocean_init_aux!(m::HBModel, p::AbstractSimpleBoxProblem, A, geom)
     A.ΔGᵘ = @SVector [-0, -0]
 
     return nothing
-
-    return nothing
 end
 
 function ocean_init_aux!(m::SWModel, p::AbstractSimpleBoxProblem, A, geom)
@@ -66,9 +60,29 @@ function ocean_init_aux!(m::SWModel, p::AbstractSimpleBoxProblem, A, geom)
     return nothing
 end
 
+"""
+    coriolis_parameter
+
+northern hemisphere coriolis
+
+# Arguments
+- `m`: model object to dispatch on and get coriolis parameters
+- `y`: y-coordinate in the box
+"""
+@inline coriolis_parameter(m::HBModel, p::AbstractSimpleBoxProblem, y) =
+    m.fₒ + m.β * y
 @inline coriolis_parameter(m::SWModel, p::AbstractSimpleBoxProblem, y) =
     m.fₒ + m.β * y
 
+############################
+# Basic box problem        #
+# Set up dimensions of box #
+############################
+
+abstract type AbstractRotation end
+struct Rotating <: AbstractRotation end
+struct Fixed <: AbstractRotation end
+
 """
     SimpleBoxProblem <: AbstractSimpleBoxProblem
 
@@ -77,7 +91,8 @@ Lˣ = zonal (east-west) length
 Lʸ = meridional (north-south) length
 H  = height of the ocean
 """
-struct SimpleBox{T, BC} <: AbstractSimpleBoxProblem
+struct SimpleBox{R, T, BC} <: AbstractSimpleBoxProblem
+    rotation::R
     Lˣ::T
     Lʸ::T
     H::T
@@ -86,60 +101,136 @@ struct SimpleBox{T, BC} <: AbstractSimpleBoxProblem
         Lˣ, # m
         Lʸ, # m
         H;  # m
+        rotation = Fixed(),
         BC = (
             OceanBC(Impenetrable(FreeSlip()), Insulating()),
             OceanBC(Penetrable(FreeSlip()), Insulating()),
         ),
     ) where {FT <: AbstractFloat}
-        return new{FT, typeof(BC)}(Lˣ, Lʸ, H, BC)
+        return new{typeof(rotation), FT, typeof(BC)}(rotation, Lˣ, Lʸ, H, BC)
     end
 end
 
-function barotropic_state!(x, t, νʰ, kˣ, gH)
-    M = @SMatrix [-νʰ * kˣ^2 gH * kˣ; -kˣ 0]
-    A = exp(M * t) * @SVector [1, 1]
+@inline coriolis_parameter(m::HBModel, ::SimpleBox{R}, y) where {R <: Fixed} =
+    -0
+@inline coriolis_parameter(m::SWModel, ::SimpleBox{R}, y) where {R <: Fixed} =
+    -0
 
-    U = A[1] * sin(kˣ * x)
-    η = A[2] * cos(kˣ * x)
+@inline coriolis_parameter(
+    m::HBModel,
+    ::SimpleBox{R},
+    y,
+) where {R <: Rotating} = m.fₒ
+@inline coriolis_parameter(
+    m::SWModel,
+    ::SimpleBox{R},
+    y,
+) where {R <: Rotating} = m.fₒ
+
+function ocean_init_state!(m::SWModel, p::SimpleBox, Q, A, coords, t)
+    k = (2π / p.Lˣ, 2π / p.Lʸ, 2π / p.H)
+    ν = (m.turbulence.ν, m.turbulence.ν, -0)
 
-    return (U = U, η = η)
+    gH = grav(m.param_set) * p.H
+    @inbounds f = coriolis_parameter(m, p, coords[2])
+
+    U, V, η = barotropic_state!(p.rotation, (coords..., t), ν, k, (gH, f))
+
+    Q.U = @SVector [U, V]
+    Q.η = η
+
+    return nothing
 end
 
 function ocean_init_state!(m::HBModel, p::SimpleBox, Q, A, coords, t)
-    @inbounds x = coords[1]
-    @inbounds y = coords[2]
-    @inbounds z = coords[3]
-
-    kˣ = 2π / p.Lˣ
-    kʸ = 2π / p.Lʸ
-    kᶻ = 2π / p.H
+    k = (2π / p.Lˣ, 2π / p.Lʸ, 2π / p.H)
+    ν = (m.νʰ, m.νʰ, m.νᶻ)
 
     gH = grav(m.param_set) * p.H
-    U, η = barotropic_state!(x, t, m.νʰ, kˣ, gH)
+    @inbounds f = coriolis_parameter(m, p, coords[2])
+
+    U, V, η = barotropic_state!(p.rotation, (coords..., t), ν, k, (gH, f))
+    u°, v° = baroclinic_deviation(p.rotation, (coords..., t), ν, k, f)
 
-    λ = m.νʰ * kˣ^2 + m.νᶻ * kᶻ^2
-    u° = exp(-λ * t) * cos(kᶻ * z) * sin(kˣ * x)
     u = u° + U / p.H
+    v = v° + V / p.H
 
-    Q.u = @SVector [u, -0]
+    Q.u = @SVector [u, v]
     Q.η = η
     Q.θ = -0
 
     return nothing
 end
 
-function ocean_init_state!(m::SWModel, p::SimpleBox, Q, A, coords, t)
-    @inbounds x = coords[1]
-    kˣ = 2π / p.Lˣ
-    νʰ = m.turbulence.ν
-    gH = grav(m.param_set) * p.H
+function barotropic_state!(
+    ::Fixed,
+    (x, y, z, t),
+    (νˣ, νʸ, νᶻ),
+    (kˣ, kʸ, kᶻ),
+    params,
+)
+    gH, _ = params
 
-    U, η = barotropic_state!(x, t, νʰ, kˣ, gH)
+    M = @SMatrix [-νˣ * kˣ^2 gH * kˣ; -kˣ 0]
+    A = exp(M * t) * @SVector [1, 1]
 
-    Q.U = @SVector [U, -0]
-    Q.η = η
+    U = A[1] * sin(kˣ * x)
+    V = -0
+    η = A[2] * cos(kˣ * x)
 
-    return nothing
+    return (U = U, V = V, η = η)
+end
+
+function baroclinic_deviation(
+    ::Fixed,
+    (x, y, z, t),
+    (νˣ, νʸ, νᶻ),
+    (kˣ, kʸ, kᶻ),
+    f,
+)
+    λ = νˣ * kˣ^2 + νᶻ * kᶻ^2
+
+    u° = exp(-λ * t) * cos(kᶻ * z) * sin(kˣ * x)
+    v° = -0
+
+    return (u° = u°, v° = v°)
+end
+
+function barotropic_state!(
+    ::Rotating,
+    (x, y, z, t),
+    (νˣ, νʸ, νᶻ),
+    (kˣ, kʸ, kᶻ),
+    params,
+)
+    gH, f = params
+
+    M = @SMatrix [-νˣ * kˣ^2 f gH * kˣ; -f -νˣ * kˣ^2 0; -kˣ 0 0]
+    A = exp(M * t) * @SVector [1, 1, 1]
+
+    U = A[1] * sin(kˣ * x)
+    V = A[2] * sin(kˣ * x)
+    η = A[3] * cos(kˣ * x)
+
+    return (U = U, V = V, η = η)
+end
+
+function baroclinic_deviation(
+    ::Rotating,
+    (x, y, z, t),
+    (νˣ, νʸ, νᶻ),
+    (kˣ, kʸ, kᶻ),
+    f,
+)
+    λ = νˣ * kˣ^2 + νᶻ * kᶻ^2
+
+    M = @SMatrix[-λ f; -f -λ]
+    A = exp(M * t) * @SVector[1, 1]
+
+    u° = A[1] * cos(kᶻ * z) * sin(kˣ * x)
+    v° = A[2] * cos(kᶻ * z) * sin(kˣ * x)
+
+    return (u° = u°, v° = v°)
 end
 
 @inline kinematic_stress(p::SimpleBox, y) = @SVector [-0, -0]
@@ -157,8 +248,6 @@ Lˣ = zonal (east-west) length
 Lʸ = meridional (north-south) length
 H  = height of the ocean
 τₒ = maximum value of wind-stress (amplitude)
-fₒ = first coriolis parameter (constant term)
-β  = second coriolis parameter (linear term)
 """
 struct HomogeneousBox{T, BC} <: AbstractSimpleBoxProblem
     Lˣ::T

src/Ocean/SplitExplicit/HydrostaticBoussinesqCoupling.jl

@@ -34,7 +34,7 @@ end
 
 @inline function coriolis_force!(m::HBModel, ::Coupled, S, Q, A, t)
     # f × u
-    f = coriolis_parameter(m, A.y)
+    f = coriolis_parameter(m, m.problem, A.y)
     uᵈ, vᵈ = A.uᵈ # Horizontal components of velocity
     S.u -= @SVector [-f * vᵈ, f * uᵈ]
 

test/Ocean/SplitExplicit/hydrostatic_spindown.jl

@@ -1,6 +1,6 @@
 include("split_explicit.jl")
 
-function SplitConfig(name, resolution, dimensions, coupling)
+function SplitConfig(name, resolution, dimensions, coupling, rotation = Fixed())
     mpicomm = MPI.COMM_WORLD
     ArrayType = ClimateMachine.array_type()
 
@@ -39,7 +39,7 @@ function SplitConfig(name, resolution, dimensions, coupling)
         polynomialorder = N,
     )
 
-    problem = SimpleBox{FT}(Lˣ, Lʸ, H)
+    problem = SimpleBox{FT}(Lˣ, Lʸ, H; rotation = rotation)
 
     model_3D = HydrostaticBoussinesqModel{FT}(
         param_set,
@@ -49,8 +49,6 @@ function SplitConfig(name, resolution, dimensions, coupling)
         αᵀ = FT(0),
         κʰ = FT(0),
         κᶻ = FT(0),
-        fₒ = FT(0),
-        β = FT(0),
     )
 
     model_2D = ShallowWaterModel{FT}(
@@ -60,8 +58,6 @@ function SplitConfig(name, resolution, dimensions, coupling)
         nothing;
         coupling = coupling,
         c = FT(1),
-        fₒ = FT(0),
-        β = FT(0),
     )
 
     return SplitConfig(

test/Ocean/SplitExplicit/test_coriolis.jl

@@ -0,0 +1,47 @@
+#!/usr/bin/env julia --project
+using Test
+
+include("hydrostatic_spindown.jl")
+ClimateMachine.init()
+
+const FT = Float64
+
+#################
+# RUN THE TESTS #
+#################
+@testset "$(@__FILE__)" begin
+
+    include("../refvals/hydrostatic_spindown_refvals.jl")
+
+    # simulation time
+    timeend = FT(15 * 24 * 3600) # s
+    tout = FT(24 * 3600) # s
+    timespan = (tout, timeend)
+
+    # DG polynomial order
+    N = Int(4)
+
+    # Domain resolution
+    Nˣ = Int(5)
+    Nʸ = Int(5)
+    Nᶻ = Int(8)
+    resolution = (N, Nˣ, Nʸ, Nᶻ)
+
+    # Domain size
+    Lˣ = 1e6  # m
+    Lʸ = 1e6  # m
+    H = 400  # m
+    dimensions = (Lˣ, Lʸ, H)
+
+    config =
+        SplitConfig("rotating", resolution, dimensions, Coupled(), Rotating())
+
+    run_split_explicit(
+        config,
+        timespan,
+        dt_fast = 300,
+        dt_slow = 300,
+        # refDat = refVals.ninety_minutes,
+        analytic_solution = true,
+    )
+end