Linear HS mountain waves (Topography)

Description of experiment

1. Dry linear Hydrostatic Mountain Waves

The atmosphere is dry and the flow impinges against a witch of Agnesi mountain of heigh $h_{m}=1$m and base parameter $a=10000m$ and centered on $x_{c} = 120km$ in a 2D domain $\Omega = 240km \times 50km$. The mountain is defined as

\[z = \frac{h_m}{1 + \frac{x - x_c}{a}}\]

The 2D problem is setup in 3D by using 1 element in the y direction. To damp the upward moving gravity waves, a Reyleigh absorbing layer is added at $z = 15000 m$.

The initial atmosphere is defined such that it has a stability frequency $N=g/\sqrt{c_p T_0}$, where

\[T_0 = 250 K\]

so that

\[\theta = \theta_0 = T_0\]

\[\pi = 1 + \frac{g^2}{c_p \theta_0 N^2}\left(\exp\left(\frac{-N^2 z}{g} \right)\right)\]

where $\theta_0 = T_0 K$.

so that

\[ρ = \frac{p_{sfc}}{R_{gas}\theta}\pi^{c_v/R_{gas}}\]

and

\[T = \theta \pi\]

1. Boundaries
   • Impenetrable(FreeSlip()) - Top and bottom: no momentum flux, no mass flux through walls.
   • Impermeable() - non-porous walls, i.e. no diffusive fluxes through walls.
   • Agnesi topography built via meshwarp.
   • Laterally periodic
2. Domain - 240,000 m (horizontal) x 4000 m (horizontal) x 30,000m (vertical)
3. Resolution - 1000m X 240 m effective resolution
4. Total simulation time - 15,000 s
5. Overrides defaults for
   • CPU Initialisation
   • Time integrator
   • Sources

Boilerplate (Using Modules)

The setup of this problem is taken from Case 6 of F.~X. { Giraldo } , M. { Restelli } (2008)

using ClimateMachine
ClimateMachine.init(parse_clargs = true);
nothing

[1638959893.639665] [hpc-92-37:26682:0]       ib_verbs.h:84   UCX  ERROR ibv_exp_query_device(mlx5_0) returned 38: No space left on device

Setting parse_clargs=true allows the use of command-line arguments (see API > Driver docs) to control simulation update and output intervals.

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using Thermodynamics.TemperatureProfiles
using Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using StaticArrays
using Test

using CLIMAParameters
using CLIMAParameters.Atmos.SubgridScale: C_smag
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

Main.##346.EarthParameterSet()

Initial Conditions

function init_agnesi_hs_lin!(problem, bl, state, aux, localgeo, t)
    (x, y, z) = localgeo.coord

    # Problem float-type
    FT = eltype(state)
    param_set = parameter_set(bl)

    # Unpack constant parameters
    R_gas::FT = R_d(param_set)
    c_p::FT = cp_d(param_set)
    c_v::FT = cv_d(param_set)
    p0::FT = MSLP(param_set)
    _grav::FT = grav(param_set)
    γ::FT = c_p / c_v

    c::FT = c_v / R_gas
    c2::FT = R_gas / c_p

    Tiso::FT = 250.0
    θ0::FT = Tiso

    # Calculate the Brunt-Vaisaila frequency for an isothermal field
    Brunt::FT = _grav / sqrt(c_p * Tiso)
    Brunt2::FT = Brunt * Brunt
    g2::FT = _grav * _grav

    π_exner::FT = exp(-_grav * z / (c_p * Tiso))
    θ::FT = θ0 * exp(Brunt2 * z / _grav)
    ρ::FT = p0 / (R_gas * θ) * (π_exner)^c

    # Compute perturbed thermodynamic state:
    T = θ * π_exner
    e_int = internal_energy(param_set, T)
    ts = PhaseDry(param_set, e_int, ρ)

    # initial velocity
    u = FT(20.0)

    # State (prognostic) variable assignment
    e_kin = FT(0)                                       # kinetic energy
    e_pot = gravitational_potential(bl.orientation, aux)# potential energy
    ρe_tot = ρ * total_energy(e_kin, e_pot, ts)         # total energy

    state.ρ = ρ
    state.ρu = SVector{3, FT}(ρ * u, 0, 0)
    state.energy.ρe = ρe_tot
end

init_agnesi_hs_lin! (generic function with 1 method)

Define a setmax method

function setmax(f, xmax, ymax, zmax)
    function setmaxima(xin, yin, zin)
        return f(xin, yin, zin; xmax = xmax, ymax = ymax, zmax = zmax)
    end
    return setmaxima
end

setmax (generic function with 1 method)

Define a warping function to build an analytic topography:

function warp_agnesi(xin, yin, zin; xmax = 1000.0, ymax = 1000.0, zmax = 1000.0)

    FT = eltype(xin)

    ac = FT(10000)
    hm = FT(1)
    xc = FT(0.5) * xmax
    zdiff = hm / (FT(1) + ((xin - xc) / ac)^2)

    # Linear relaxation towards domain maximum height
    x, y, z = xin, yin, zin + zdiff * (zmax - zin) / zmax
    return x, y, z
end

warp_agnesi (generic function with 1 method)

Model Configuration

We define a configuration function to assist in prescribing the physical model. The purpose of this is to populate the AtmosLESConfiguration with arguments appropriate to the problem being considered.

function config_agnesi_hs_lin(
    ::Type{FT},
    N,
    resolution,
    xmax,
    ymax,
    zmax,
) where {FT}
    #
    # Explicit Rayleigh damping:
    #
    # ``
    #   \tau_s = \alpha * \sin\left(0.5\pi \frac{z - z_s}{zmax - z_s} \right)^2,
    # ``
    # where
    # ``sponge_ampz`` is the wave damping coefficient (1/s)
    # ``z_s`` is the level where the Rayleigh sponge starts
    # ``zmax`` is the domain top
    #
    # Setup the parameters for the gravity wave absorbing layer
    # at the top of the domain
    #
    # u_relaxation(xvelo, vvelo, wvelo) contains the background velocity values to which
    # the sponge relaxes the vertically moving wave
    u_relaxation = SVector(FT(20), FT(0), FT(0))

    # Wave damping coefficient (1/s)
    sponge_ampz = FT(0.5)

    # Vertical level where the absorbing layer starts
    z_s = FT(25000.0)

    # Pass the sponge parameters to the sponge calculator
    rayleigh_sponge =
        RayleighSponge{FT}(zmax, z_s, sponge_ampz, u_relaxation, 2)

    # Setup the source terms for this problem:
    source = (Gravity(), rayleigh_sponge)

    # Define the reference state:
    T_virt = FT(250)
    temp_profile_ref = IsothermalProfile(param_set, T_virt)
    ref_state = HydrostaticState(temp_profile_ref)
    nothing # hide

    _C_smag = FT(0.21)
    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = Vreman(_C_smag),
        moisture = DryModel(),
        tracers = NoTracers(),
    )
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        physics;
        init_state_prognostic = init_agnesi_hs_lin!,
        source = source,
    )

    config = ClimateMachine.AtmosLESConfiguration(
        "Agnesi_HS_LINEAR",      # Problem title [String]
        N,                       # Polynomial order [Int]
        resolution,              # (Δx, Δy, Δz) effective resolution [m]
        xmax,                    # Domain maximum size [m]
        ymax,                    # Domain maximum size [m]
        zmax,                    # Domain maximum size [m]
        param_set,               # Parameter set.
        init_agnesi_hs_lin!,     # Function specifying initial condition
        model = model,           # Model type
        meshwarp = setmax(warp_agnesi, xmax, ymax, zmax),
    )

    return config
end

config_agnesi_hs_lin (generic function with 1 method)

Define a main method (entry point)

function main()

    FT = Float64

    # Define the polynomial order and effective grid spacings:
    N = 4

    # Define the domain size and spatial resolution
    Nx = 20
    Ny = 20
    Nz = 20
    xmax = FT(244000)
    ymax = FT(4000)
    zmax = FT(50000)
    Δx = xmax / FT(Nx)
    Δy = ymax / FT(Ny)
    Δz = zmax / FT(Nz)
    resolution = (Δx, Δy, Δz)

    t0 = FT(0)
    timeend = FT(150) #FT(hrs * 60 * 60)

    # Define the max Courant for the time time integrator (ode_solver).
    # The default value is 1.7 for LSRK144:
    CFL = FT(1.5)

    # Assign configurations so they can be passed to the `invoke!` function
    driver_config = config_agnesi_hs_lin(FT, N, resolution, xmax, ymax, zmax)

    # Define the time integrator:
    # We chose an explicit single-rate LSRK144 for this problem
    ode_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK144NiegemannDiehlBusch,
    )

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        init_on_cpu = true,
        Courant_number = CFL,
    )

    # Set up the spectral filter to remove the solutions spurious modes
    # Define the order of the exponential filter: use 32 or 64 for this problem.
    # The larger the value, the less dissipation you get:
    filterorder = 64
    filter = ExponentialFilter(solver_config.dg.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            AtmosFilterPerturbations(driver_config.bl),
            solver_config.dg.grid,
            filter,
            state_auxiliary = solver_config.dg.state_auxiliary,
        )
        nothing
    end
    # End exponential filter

    # Invoke solver (calls `solve!` function for time-integrator),
    # pass the driver, solver and diagnostic config information.
    result = ClimateMachine.invoke!(
        solver_config;
        user_callbacks = (cbfilter,),
        check_euclidean_distance = true,
    )

    # Check that the solution norm is reasonable.
    @test isapprox(result, FT(1); atol = 1.5e-3)
end

main (generic function with 1 method)

Call main

main()

Test Passed

This page was generated using Literate.jl.