From 94c53eb2e836f45e8ff35b6b0101a0f5cd084ef5 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Thu, 11 Jul 2024 11:19:50 -0800 Subject: [PATCH] Improve notation for computing turbulent fluxes (#99) * Improve notation for computing turbulent fluxes * Minor cleanup * Add some comments * Another comment * Update src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl Co-authored-by: Simone Silvestri * Resolve dependencies * Add using ClimaOcean to runtests setup * Fix bug in extending - --------- Co-authored-by: Simone Silvestri --- Manifest.toml | 6 +- .../atmosphere_ocean_fluxes.jl | 60 +++++++------ .../ocean_sea_ice_surface_fluxes.jl | 2 + .../similarity_theory_turbulent_fluxes.jl | 87 ++++++++++--------- .../CrossRealmFluxes/stability_functions.jl | 15 ++-- test/runtests_setup.jl | 2 + 6 files changed, 93 insertions(+), 79 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index aab4001a..d88bbd31 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.4" +julia_version = "1.10.0" manifest_format = "2.0" project_hash = "086b5c52c3f2a61f133762647db63128f8714c56" @@ -255,7 +255,7 @@ weakdeps = ["Dates", "LinearAlgebra"] [[deps.CompilerSupportLibraries_jll]] deps = ["Artifacts", "Libdl"] uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae" -version = "1.1.1+0" +version = "1.0.5+1" [[deps.CompositionsBase]] git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad" @@ -906,7 +906,7 @@ weakdeps = ["Adapt"] [[deps.OpenBLAS_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"] uuid = "4536629a-c528-5b80-bd46-f80d51c5b363" -version = "0.3.23+4" +version = "0.3.23+2" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl index 7c1c9f44..2c9c368f 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/atmosphere_ocean_fluxes.jl @@ -58,7 +58,8 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) # kernel parameters that compute fluxes in 0:Nx+1 and 0:Ny+1 kernel_parameters = KernelParameters(kernel_size, (-1, -1)) - launch!(arch, grid, kernel_parameters, _compute_atmosphere_ocean_similarity_theory_fluxes!, + launch!(arch, grid, kernel_parameters, + _compute_atmosphere_ocean_similarity_theory_fluxes!, similarity_theory, grid, clock, @@ -73,7 +74,8 @@ function compute_atmosphere_ocean_fluxes!(coupled_model) atmosphere.boundary_layer_height, atmosphere.thermodynamics_parameters) - launch!(arch, grid, kernel_parameters, _assemble_atmosphere_ocean_fluxes!, + launch!(arch, grid, kernel_parameters, + _assemble_atmosphere_ocean_fluxes!, centered_velocity_fluxes, net_tracer_fluxes, grid, @@ -182,19 +184,13 @@ limit_fluxes_over_sea_ice!(args...) = nothing similarity_theory.water_vapor_saturation, surface_type) - # Thermodynamic and dynamic surface state + # Thermodynamic and dynamic (ocean) surface state 𝒬₀ = thermodynamic_surface_state = AtmosphericThermodynamics.PhaseEquil_pTq(ℂₐ, pₐ, Tₒ, qₒ) - h₀ = zero(grid) # surface height Uₒ = SVector(uₒ, vₒ) 𝒰₀ = dynamic_ocean_state = SurfaceFluxes.StateValues(h₀, Uₒ, 𝒬₀) - Qv = similarity_theory.fields.latent_heat - Qc = similarity_theory.fields.sensible_heat - Fv = similarity_theory.fields.water_vapor - τx = similarity_theory.fields.x_momentum - τy = similarity_theory.fields.y_momentum - + # Some parameters g = default_gravitational_acceleration ϰ = similarity_theory.von_karman_constant @@ -209,16 +205,24 @@ limit_fluxes_over_sea_ice!(args...) = nothing # Convert back from a zonal - meridional flux to the frame of # reference of the native ocean grid - τˣ, τʸ = convert_to_extrinsic_reference_frame(i, j, kᴺ, grid, turbulent_fluxes.x_momentum, - turbulent_fluxes.y_momentum) + ρτxⁱʲ, ρτyⁱʲ = convert_to_extrinsic_reference_frame(i, j, kᴺ, grid, + turbulent_fluxes.x_momentum, + turbulent_fluxes.y_momentum) + + # Store fluxes + Qv = similarity_theory.fields.latent_heat + Qc = similarity_theory.fields.sensible_heat + Fv = similarity_theory.fields.water_vapor + ρτx = similarity_theory.fields.x_momentum + ρτy = similarity_theory.fields.y_momentum @inbounds begin # +0: cooling, -0: heating - Qv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.latent_heat) - Qc[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.sensible_heat) - Fv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.water_vapor) - τx[i, j, 1] = ifelse(inactive, 0, τˣ) - τy[i, j, 1] = ifelse(inactive, 0, τʸ) + Qv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.latent_heat) + Qc[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.sensible_heat) + Fv[i, j, 1] = ifelse(inactive, 0, turbulent_fluxes.water_vapor) + ρτx[i, j, 1] = ifelse(inactive, 0, ρτxⁱʲ) + ρτy[i, j, 1] = ifelse(inactive, 0, ρτyⁱʲ) end end @@ -262,11 +266,11 @@ end Mp = interp_atmos_time_series(prescribed_freshwater_flux, X, time, atmos_args...) Mr = get_runoff_flux(X, time, runoff_args) - Qc = similarity_theory_fields.sensible_heat[i, j, 1] # sensible or "conductive" heat flux - Qv = similarity_theory_fields.latent_heat[i, j, 1] # latent heat flux - Mv = similarity_theory_fields.water_vapor[i, j, 1] # mass flux of water vapor - τx = similarity_theory_fields.x_momentum[i, j, 1] # zonal momentum flux - τy = similarity_theory_fields.y_momentum[i, j, 1] # meridional momentum flux + Qc = similarity_theory_fields.sensible_heat[i, j, 1] # sensible or "conductive" heat flux + Qv = similarity_theory_fields.latent_heat[i, j, 1] # latent heat flux + Mv = similarity_theory_fields.water_vapor[i, j, 1] # mass flux of water vapor + ρτx = similarity_theory_fields.x_momentum[i, j, 1] # zonal momentum flux + ρτy = similarity_theory_fields.y_momentum[i, j, 1] # meridional momentum flux end # Compute heat fluxes, bulk flux first @@ -286,16 +290,16 @@ end ΣF += Fv # Compute fluxes for u, v, T, S from momentum, heat, and freshwater fluxes - Jᵘ = centered_velocity_fluxes.u - Jᵛ = centered_velocity_fluxes.v + τx = centered_velocity_fluxes.u + τy = centered_velocity_fluxes.v Jᵀ = net_tracer_fluxes.T Jˢ = net_tracer_fluxes.S ρₒ = ocean_reference_density cₒ = ocean_heat_capacity - atmos_ocean_Jᵘ = τx / ρₒ - atmos_ocean_Jᵛ = τy / ρₒ + atmos_ocean_τx = ρτx / ρₒ + atmos_ocean_τy = ρτy / ρₒ atmos_ocean_Jᵀ = ΣQ / (ρₒ * cₒ) atmos_ocean_Jˢ = - Sₒ * ΣF @@ -303,8 +307,8 @@ end inactive = inactive_node(i, j, kᴺ, grid, c, c, c) @inbounds begin - Jᵘ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵘ) - Jᵛ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵛ) + τx[i, j, 1] = ifelse(inactive, 0, atmos_ocean_τx) + τy[i, j, 1] = ifelse(inactive, 0, atmos_ocean_τy) Jᵀ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jᵀ) Jˢ[i, j, 1] = ifelse(inactive, 0, atmos_ocean_Jˢ) end diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/ocean_sea_ice_surface_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/ocean_sea_ice_surface_fluxes.jl index 61899ac1..6c8dcb81 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/ocean_sea_ice_surface_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/ocean_sea_ice_surface_fluxes.jl @@ -29,6 +29,8 @@ using KernelAbstractions: @kernel, @index struct OceanSeaIceSurfaceFluxes{T, P, C, R, PI, PC, FT, UN} turbulent :: T prescribed :: P + # Add `components` which will also store components of the total fluxes + # (eg latent, sensible heat flux) total :: C radiation :: R previous_ice_thickness :: PI diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl index 7cb5ebf3..c1151a8f 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/similarity_theory_turbulent_fluxes.jl @@ -30,20 +30,20 @@ import SurfaceFluxes.Parameters: ##### struct SimilarityTheoryTurbulentFluxes{FT, UF, TP, S, W, R, B, V, F} - gravitational_acceleration :: FT - von_karman_constant :: FT - turbulent_prandtl_number :: FT - gustiness_parameter :: FT - stability_functions :: UF - thermodynamics_parameters :: TP - water_vapor_saturation :: S - water_mole_fraction :: W - roughness_lengths :: R - bulk_coefficients :: B - bulk_velocity :: V - tolerance :: FT - maxiter :: Int - fields :: F + gravitational_acceleration :: FT # parameter + von_karman_constant :: FT # parameter + turbulent_prandtl_number :: FT # parameter + gustiness_parameter :: FT # bulk velocity parameter + stability_functions :: UF # functions for turbulent fluxes + thermodynamics_parameters :: TP # parameter group + water_vapor_saturation :: S # model for computing the saturation water vapor mass + water_mole_fraction :: W # mole fraction of H₂O in seawater + roughness_lengths :: R # parameterization for turbulent fluxes + bulk_coefficients :: B # ? + bulk_velocity :: V # bulk velocity scale for turbulent fluxes + tolerance :: FT # solver option + maxiter :: Int # solver option + fields :: F # fields that store turbulent fluxes end const STTF = SimilarityTheoryTurbulentFluxes @@ -229,23 +229,27 @@ end Σ★ = SimilarityScales(u★, u★, u★) # The inital velocity scale assumes that - # the gustiness velocity `uᴳ` is equal to 0.5 ms⁻¹. + # the gustiness velocity `Uᴳ` is equal to 0.5 ms⁻¹. # That will be refined later on. - ΔUᴳ = sqrt(Δu^2 + Δv^2 + convert(eltype(Δh), 0.25)) + FT = eltype(Δh) + Uᴳᵢ = convert(FT, 0.5^2) + ΔU = sqrt(Δu^2 + Δv^2 + Uᴳᵢ) # Initialize the solver iteration = 0 while iterating(Σ★ - Σ₀, iteration, maxiter, similarity_theory) Σ₀ = Σ★ - Σ★, ΔUᴳ = refine_characteristic_scales(Σ★, ΔUᴳ, - similarity_theory, - surface_state, - differences, - atmos_boundary_layer_height, - thermodynamics_parameters, - gravitational_acceleration, - von_karman_constant) + # Refine both the characteristic scale and the effective + # velocity difference ΔU, including gustiness. + Σ★, ΔU = refine_similarity_variables(Σ★, ΔU, + similarity_theory, + surface_state, + differences, + atmos_boundary_layer_height, + thermodynamics_parameters, + gravitational_acceleration, + von_karman_constant) iteration += 1 end @@ -257,9 +261,9 @@ end q★ = q★ / similarity_theory.turbulent_prandtl_number # `u★² ≡ sqrt(τx² + τy²)` - # We remove the gustiness by dividing by `ΔUᴳ` - τx = - u★^2 * Δu / ΔUᴳ - τy = - u★^2 * Δv / ΔUᴳ + # We remove the gustiness by dividing by `ΔU` + τx = - u★^2 * Δu / ΔU + τy = - u★^2 * Δv / ΔU 𝒬ₐ = atmos_state.ts ρₐ = AtmosphericThermodynamics.air_density(ℂₐ, 𝒬ₐ) @@ -326,15 +330,15 @@ end return Δh, Δu, Δv, Δθ, Δq end -@inline function refine_characteristic_scales(estimated_characteristic_scales, - velocity_scale, - similarity_theory, - surface_state, - differences, - atmos_boundary_layer_height, - thermodynamics_parameters, - gravitational_acceleration, - von_karman_constant) +@inline function refine_similarity_variables(estimated_characteristic_scales, + velocity_scale, + similarity_theory, + surface_state, + differences, + atmos_boundary_layer_height, + thermodynamics_parameters, + gravitational_acceleration, + von_karman_constant) # "initial" scales because we will recompute them u★ = estimated_characteristic_scales.momentum @@ -358,7 +362,7 @@ end ℂ = thermodynamics_parameters g = gravitational_acceleration 𝒬ₒ = surface_state.ts # thermodynamic state - zᵢ = atmos_boundary_layer_height + hᵢ = atmos_boundary_layer_height # Compute Monin-Obukhov length scale depending on a `buoyancy flux` b★ = buoyancy_scale(θ★, q★, 𝒬ₒ, ℂ, g) @@ -388,10 +392,11 @@ end # Buoyancy flux characteristic scale for gustiness (Edson 2013) Jᵇ = - u★ * b★ - uᴳ = β * cbrt(Jᵇ * zᵢ) + Uᴳ = β * cbrt(Jᵇ * hᵢ) # New velocity difference accounting for gustiness - ΔUᴳ = sqrt(Δu^2 + Δv^2 + uᴳ^2) + ΔU = sqrt(Δu^2 + Δv^2 + Uᴳ^2) + + return SimilarityScales(u★, θ★, q★), ΔU +end - return SimilarityScales(u★, θ★, q★), ΔUᴳ -end \ No newline at end of file diff --git a/src/OceanSeaIceModels/CrossRealmFluxes/stability_functions.jl b/src/OceanSeaIceModels/CrossRealmFluxes/stability_functions.jl index f1f67cdf..b30845e2 100644 --- a/src/OceanSeaIceModels/CrossRealmFluxes/stability_functions.jl +++ b/src/OceanSeaIceModels/CrossRealmFluxes/stability_functions.jl @@ -1,5 +1,5 @@ -import Base: - import Statistics +import Base: - ##### ##### Struct that represents a 3-tuple of momentum, heat, and water vapor @@ -11,18 +11,19 @@ struct SimilarityScales{U, T, Q} water_vapor :: Q end --(a::SimilarityScales, b::SimilarityScales) = SimilarityScales(a.momentum - b.momentum, - a.temperature - b.temperature, - a.water_vapor - b.water_vapor) +function -(a::SimilarityScales, b::SimilarityScales) + Δu = a.momentum - b.momentum + Δθ = a.temperature - b.temperature + Δq = a.water_vapor - b.water_vapor + return SimilarityScales(Δu, Δθ, Δq) +end Statistics.norm(a::SimilarityScales) = norm(a.momentum) + norm(a.temperature) + norm(a.water_vapor) +# Edson et al. (2013) function edson_stability_functions(FT = Float64) - - # Edson et al. (2013) ψu = MomentumStabilityFunction() ψc = ScalarStabilityFunction() - return SimilarityScales(ψu, ψc, ψc) end diff --git a/test/runtests_setup.jl b/test/runtests_setup.jl index 83e3f1d9..d1cca936 100644 --- a/test/runtests_setup.jl +++ b/test/runtests_setup.jl @@ -6,4 +6,6 @@ using Test using Oceananigans.Architectures: architecture using Oceananigans.OutputReaders: interpolate! +using ClimaOcean + test_architectures = [CPU()]