From 170cf67624bc52415970feac4881fe629a0cb034 Mon Sep 17 00:00:00 2001 From: yairchn Date: Sun, 3 Oct 2021 08:55:53 -0700 Subject: [PATCH] explict gmv equations --- integration_tests/utils/generate_namelist.jl | 8 +- src/Operators.jl | 20 +++ src/Turbulence_PrognosticTKE.jl | 124 +++++++++---------- src/Variables.jl | 11 -- src/types.jl | 12 +- 5 files changed, 87 insertions(+), 88 deletions(-) diff --git a/integration_tests/utils/generate_namelist.jl b/integration_tests/utils/generate_namelist.jl index 67802b9262..d174dca904 100644 --- a/integration_tests/utils/generate_namelist.jl +++ b/integration_tests/utils/generate_namelist.jl @@ -128,7 +128,7 @@ function default_namelist(case_name::String) namelist_defaults["thermodynamics"]["quadrature_type"] = "log-normal" #"gaussian" or "log-normal" namelist_defaults["time_stepping"] = Dict() - namelist_defaults["time_stepping"]["dt"] = 10.0 + namelist_defaults["time_stepping"]["dt"] = 7.0 namelist_defaults["microphysics"] = Dict() namelist_defaults["microphysics"]["rain_model"] = "None" @@ -351,11 +351,11 @@ function SP(namelist_defaults) namelist["meta"]["casename"] = "SP" # this case is resolution dependent, we should check why - namelist["grid"]["nz"] = 256 - namelist["grid"]["dz"] = 8 + namelist["grid"]["nz"] = 128 + namelist["grid"]["dz"] = 16 namelist["time_stepping"]["t_max"] = 7200.0 - namelist["time_stepping"]["dt"] = 3.0 + namelist["time_stepping"]["dt"] = 1.0 namelist["meta"]["simname"] = "SP" namelist["meta"]["casename"] = "SP" diff --git a/src/Operators.jl b/src/Operators.jl index 52ac912e16..bb9d263031 100644 --- a/src/Operators.jl +++ b/src/Operators.jl @@ -36,6 +36,26 @@ end ∇f2c(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[2] - bc.value) * grid.Δzi ∇f2c(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value +∇c2f(f_dual::SA.SVector, grid::Grid, k::Int; bottom = UseBoundaryValue(), top = UseBoundaryValue()) = + ∇c2f(f_dual, grid, k, bottom, top) + +function ∇c2f(f_dual::SA.SVector, grid::Grid, k::Int, bottom::AbstractBC, top::AbstractBC) + if is_surface_face(grid, k) + return ∇c2f(f_dual, grid, BottomBCTag(), bottom) + elseif is_toa_face(grid, k) + return ∇c2f(f_dual, grid, TopBCTag(), top) + else + return ∇c2f(f_dual, grid, InteriorTag()) + end +end +∇c2f(f::SA.SVector, grid::Grid, ::Int, ::UseBoundaryValue, top::UseBoundaryValue) = ∇c2f(f, grid, InteriorTag()) +∇c2f(f::SA.SVector, grid::Grid, ::InteriorTag) = (f[2] - f[1]) * grid.Δzi +∇c2f(f::SA.SVector, grid::Grid, ::TopBCTag, bc::SetValue) = (bc.value - f[1]) * grid.Δzi * 2.0 +∇c2f(f::SA.SVector, grid::Grid, ::TopBCTag, bc::SetGradient) = bc.value +∇c2f(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetValue) = (f[2] - bc.value) * grid.Δzi * 2.0 +∇c2f(f::SA.SVector, grid::Grid, ::BottomBCTag, bc::SetGradient) = bc.value + + # TODO: this should be changed to `c∇` or `c∇_upwind` # Actually, this method is needed for rain (upwind is in reverse direction due to rain) function c∇_downwind(f_dual::SA.SVector, grid::Grid, k::Int; bottom = NoBCGivenError(), top = NoBCGivenError()) diff --git a/src/Turbulence_PrognosticTKE.jl b/src/Turbulence_PrognosticTKE.jl index 1a4d8c4403..313c66b912 100755 --- a/src/Turbulence_PrognosticTKE.jl +++ b/src/Turbulence_PrognosticTKE.jl @@ -51,10 +51,10 @@ function initialize_io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats) add_profile(Stats, "massflux_qt") add_profile(Stats, "massflux_tendency_h") add_profile(Stats, "massflux_tendency_qt") - add_profile(Stats, "diffusive_flux_h") - add_profile(Stats, "diffusive_flux_u") - add_profile(Stats, "diffusive_flux_v") - add_profile(Stats, "diffusive_flux_qt") + # add_profile(Stats, "diffusive_flux_h") + # add_profile(Stats, "diffusive_flux_u") + # add_profile(Stats, "diffusive_flux_v") + # add_profile(Stats, "diffusive_flux_qt") add_profile(Stats, "diffusive_tendency_h") add_profile(Stats, "diffusive_tendency_qt") add_profile(Stats, "total_flux_h") @@ -173,14 +173,14 @@ function io(edmf::EDMF_PrognosticTKE, Stats::NetCDFIO_Stats, TS::TimeStepping) write_profile(Stats, "massflux_qt", mf_qt) write_profile(Stats, "massflux_tendency_h", edmf.massflux_tendency_h) write_profile(Stats, "massflux_tendency_qt", edmf.massflux_tendency_qt) - write_profile(Stats, "diffusive_flux_h", edmf.diffusive_flux_h) - write_profile(Stats, "diffusive_flux_qt", edmf.diffusive_flux_qt) - write_profile(Stats, "diffusive_flux_u", edmf.diffusive_flux_u) - write_profile(Stats, "diffusive_flux_v", edmf.diffusive_flux_v) + # write_profile(Stats, "diffusive_flux_h", edmf.diffusive_flux_h) + # write_profile(Stats, "diffusive_flux_qt", edmf.diffusive_flux_qt) + # write_profile(Stats, "diffusive_flux_u", edmf.diffusive_flux_u) + # write_profile(Stats, "diffusive_flux_v", edmf.diffusive_flux_v) write_profile(Stats, "diffusive_tendency_h", edmf.diffusive_tendency_h) write_profile(Stats, "diffusive_tendency_qt", edmf.diffusive_tendency_qt) - write_profile(Stats, "total_flux_h", mf_h .+ edmf.diffusive_flux_h) - write_profile(Stats, "total_flux_qt", mf_qt .+ edmf.diffusive_flux_qt) + # write_profile(Stats, "total_flux_h", mf_h .+ interpf2c.(edmf.diffusive_flux_h)) + # write_profile(Stats, "total_flux_qt", mf_qt .+ interpf2c.(edmf.diffusive_flux_qt)) write_profile(Stats, "mixing_length", edmf.mixing_length) write_profile(Stats, "updraft_qt_precip", edmf.UpdThermo.prec_source_qt_tot) write_profile(Stats, "updraft_thetal_precip", edmf.UpdThermo.prec_source_h_tot) @@ -254,6 +254,7 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, Case, GMV, ref_s param_set = parameter_set(GMV) rho0 = ref_state.rho0 kf_surf = kf_surface(grid) + kc_surf = kc_surface(grid) ae = 1 .- edmf.UpdVar.Area.bulkvalues # area of environment @inbounds for k in real_center_indices(grid) @@ -369,6 +370,28 @@ function compute_gm_tendencies!(edmf::EDMF_PrognosticTKE, grid, Case, GMV, ref_s GMV.H.tendencies[k] += mf_tend_h GMV.QT.tendencies[k] += mf_tend_qt end + aeKHq_tot_bc = -Case.Sur.rho_qtflux / edmf.ae[kc_surf] + aeKHθ_liq_ice_bc = -Case.Sur.rho_hflux / edmf.ae[kc_surf] + aeKHu_bc = -Case.Sur.rho_uflux / edmf.ae[kc_surf] + aeKHv_bc = -Case.Sur.rho_vflux / edmf.ae[kc_surf] + @inbounds for k in real_center_indices(grid) + α0_half = ref_state.alpha0_half[k] + aeKH_q_tot_cut = dual_faces(edmf.diffusive_flux_qt, grid, k) + ∇aeKH_q_tot = ∇f2c(aeKH_q_tot_cut, grid, k; bottom = SetValue(aeKHq_tot_bc), top = SetValue(0)) + GMV.QT.tendencies[k] += α0_half * edmf.ae[k] * ∇aeKH_q_tot + + aeKH_θ_liq_ice_cut = dual_faces(edmf.diffusive_flux_h, grid, k) + ∇aeKH_θ_liq_ice = ∇f2c(aeKH_θ_liq_ice_cut, grid, k; bottom = SetValue(aeKHθ_liq_ice_bc), top = SetValue(0)) + GMV.H.tendencies[k] += α0_half * edmf.ae[k] * ∇aeKH_θ_liq_ice + + aeKM_u_cut = dual_faces(edmf.diffusive_flux_u, grid, k) + ∇aeKM_u = ∇f2c(aeKM_u_cut, grid, k; bottom = SetValue(aeKHu_bc), top = SetValue(0)) + GMV.U.tendencies[k] += α0_half * edmf.ae[k] * ∇aeKM_u + + aeKM_v_cut = dual_faces(edmf.diffusive_flux_v, grid, k) + ∇aeKM_v = ∇f2c(aeKM_v_cut, grid, k; bottom = SetValue(aeKHv_bc), top = SetValue(0)) + GMV.V.tendencies[k] += α0_half * edmf.ae[k] * ∇aeKM_v + end end function compute_diffusive_fluxes(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase, TS::TimeStepping) @@ -390,26 +413,26 @@ function compute_diffusive_fluxes(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariabl edmf.rho_ae_KH[k] = interpc2f(aeKH, grid, k; aeKH_bcs...) * ref_state.rho0[k] edmf.rho_ae_KM[k] = interpc2f(aeKM, grid, k; aeKM_bcs...) * ref_state.rho0[k] end - q_bc = -Case.Sur.rho_qtflux / edmf.rho_ae_KH[kf_surf] - θ_liq_ice_bc = -Case.Sur.rho_hflux / edmf.rho_ae_KH[kf_surf] - u_bc = -Case.Sur.rho_uflux / edmf.rho_ae_KM[kf_surf] - v_bc = -Case.Sur.rho_vflux / edmf.rho_ae_KM[kf_surf] - @inbounds for k in real_center_indices(grid) - q_cut = ccut(edmf.EnvVar.QT.values, grid, k) - ∇q_tot = c∇(q_cut, grid, k; bottom = SetGradient(q_bc), top = SetGradient(0)) - edmf.diffusive_flux_qt[k] = -0.5 * ref_state.rho0_half[k] * edmf.ae[k] * KH[k] * ∇q_tot + q_bc = Case.Sur.rho_qtflux / edmf.ae[kc_surf] + θ_liq_ice_bc = Case.Sur.rho_hflux / edmf.ae[kc_surf] + u_bc = Case.Sur.rho_uflux / edmf.ae[kc_surf] + v_bc = Case.Sur.rho_vflux / edmf.ae[kc_surf] + @inbounds for k in real_face_indices(grid) + q_dual = dual_centers(edmf.EnvVar.QT.values, grid, k) + ∇q_tot_face = ∇c2f(q_dual, grid, k; bottom = SetGradient(q_bc), top = SetGradient(0)) + edmf.diffusive_flux_qt[k] = edmf.rho_ae_KH[k] * ∇q_tot_face - θ_liq_ice_cut = ccut(edmf.EnvVar.H.values, grid, k) - ∇θ_liq_ice = c∇(θ_liq_ice_cut, grid, k; bottom = SetGradient(θ_liq_ice_bc), top = SetGradient(0)) - edmf.diffusive_flux_h[k] = -0.5 * ref_state.rho0_half[k] * edmf.ae[k] * KH[k] * ∇θ_liq_ice + θ_liq_ice_dual = dual_centers(edmf.EnvVar.H.values, grid, k) + ∇θ_liq_ice_face = ∇c2f(θ_liq_ice_dual, grid, k; bottom = SetGradient(θ_liq_ice_bc), top = SetGradient(0)) + edmf.diffusive_flux_h[k] = edmf.rho_ae_KH[k] * ∇θ_liq_ice_face - u_cut = ccut(GMV.U.values, grid, k) - ∇u = c∇(u_cut, grid, k; bottom = SetGradient(u_bc), top = SetGradient(0)) - edmf.diffusive_flux_u[k] = -0.5 * ref_state.rho0_half[k] * edmf.ae[k] * KM[k] * ∇u + u_face = dual_centers(GMV.U.values, grid, k) + ∇u_face = ∇c2f(u_face, grid, k; bottom = SetGradient(u_bc), top = SetGradient(0)) + edmf.diffusive_flux_u[k] = edmf.rho_ae_KM[k] * ∇u_face - v_cut = ccut(GMV.V.values, grid, k) - ∇v = c∇(v_cut, grid, k; bottom = SetGradient(v_bc), top = SetGradient(0)) - edmf.diffusive_flux_v[k] = -0.5 * ref_state.rho0_half[k] * edmf.ae[k] * KM[k] * ∇v + v_face = dual_centers(GMV.V.values, grid, k) + ∇v_face = ∇c2f(v_face, grid, k; bottom = SetGradient(v_bc), top = SetGradient(0)) + edmf.diffusive_flux_v[k] = edmf.rho_ae_KM[k] * ∇v_face end return end @@ -592,21 +615,6 @@ function update(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBas # ----------- TODO: move to compute_tendencies implicit_eqs = edmf.implicit_eqs # Matrix is the same for all variables that use the same eddy diffusivity, we can construct once and reuse - implicit_eqs.A_θq_gm .= construct_tridiag_diffusion_gm(grid, TS.dt, edmf.rho_ae_KH, ref_state.rho0_half, edmf.ae) - implicit_eqs.A_uv_gm .= construct_tridiag_diffusion_gm(grid, TS.dt, edmf.rho_ae_KM, ref_state.rho0_half, edmf.ae) - Δzi = grid.Δzi - @inbounds for k in real_center_indices(grid) - implicit_eqs.b_u_gm[k] = GMV.U.values[k] - implicit_eqs.b_v_gm[k] = GMV.V.values[k] - implicit_eqs.b_q_tot_gm[k] = edmf.EnvVar.QT.values[k] - implicit_eqs.b_θ_liq_ice_gm[k] = edmf.EnvVar.H.values[k] - if is_surface_center(grid, k) - implicit_eqs.b_u_gm[k] += TS.dt * Case.Sur.rho_uflux * Δzi * ref_state.alpha0_half[k] / edmf.ae[k] - implicit_eqs.b_v_gm[k] += TS.dt * Case.Sur.rho_vflux * Δzi * ref_state.alpha0_half[k] / edmf.ae[k] - implicit_eqs.b_q_tot_gm[k] += TS.dt * Case.Sur.rho_qtflux * Δzi * ref_state.alpha0_half[k] / edmf.ae[k] - implicit_eqs.b_θ_liq_ice_gm[k] += TS.dt * Case.Sur.rho_hflux * Δzi * ref_state.alpha0_half[k] / edmf.ae[k] - end - end gm = GMV up = edmf.UpdVar @@ -646,25 +654,11 @@ function update(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBas # Update the grid mean variables with the tendency due to eddy diffusion # Solve tri-diagonal systems - x_q_tot_gm = center_field(grid) # for tridiag solver - x_θ_liq_ice_gm = center_field(grid) # for tridiag solver - x_q_tot_gm .= tridiag_solve(implicit_eqs.b_q_tot_gm, implicit_eqs.A_θq_gm) - x_θ_liq_ice_gm .= tridiag_solve(implicit_eqs.b_θ_liq_ice_gm, implicit_eqs.A_θq_gm) - GMV.U.new .= tridiag_solve(implicit_eqs.b_u_gm, implicit_eqs.A_uv_gm) - GMV.V.new .= tridiag_solve(implicit_eqs.b_v_gm, implicit_eqs.A_uv_gm) en.TKE.values .= tridiag_solve(implicit_eqs.b_TKE, implicit_eqs.A_TKE) en.Hvar.values .= tridiag_solve(implicit_eqs.b_Hvar, implicit_eqs.A_Hvar) en.QTvar.values .= tridiag_solve(implicit_eqs.b_QTvar, implicit_eqs.A_QTvar) en.HQTcov.values .= tridiag_solve(implicit_eqs.b_HQTcov, implicit_eqs.A_HQTcov) - @inbounds for k in real_center_indices(grid) - GMV.QT.new[k] = max(GMV.QT.values[k] + edmf.ae[k] * (x_q_tot_gm[k] - edmf.EnvVar.QT.values[k]), 0.0) - GMV.H.new[k] = GMV.H.values[k] + edmf.ae[k] * (x_θ_liq_ice_gm[k] - edmf.EnvVar.H.values[k]) - # get the diffusive flux, TODO: move to diagnostics (callbacks?) - edmf.diffusive_tendency_h[k] = (GMV.H.new[k] - GMV.H.values[k]) * TS.dti - edmf.diffusive_tendency_qt[k] = (GMV.QT.new[k] - GMV.QT.values[k]) * TS.dti - end - # Filter solution, TODO: fuse with `zero_area_fraction_cleanup` and put into `filter_variables!` @inbounds for k in real_center_indices(grid) en.TKE.values[k] = max(en.TKE.values[k], 0.0) @@ -674,27 +668,25 @@ function update(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBas en.HQTcov.values[k] = min(en.HQTcov.values[k], sqrt(en.Hvar.values[k] * en.QTvar.values[k])) end - update_GMV_turbulence(edmf, GMV, Case, TS) # set values set_values_with_new(edmf.UpdVar) zero_area_fraction_cleanup(edmf, GMV) - update(GMV, TS) + GMV_update(grid, GMV, TS) return end -function update_GMV_turbulence(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase, TS::TimeStepping) - - @inbounds for k in real_center_indices(edmf.grid) - GMV.H.tendencies[k] += (GMV.H.new[k] - GMV.H.values[k]) * TS.dti - GMV.QT.tendencies[k] += (GMV.QT.new[k] - GMV.QT.values[k]) * TS.dti - GMV.U.tendencies[k] += (GMV.U.new[k] - GMV.U.values[k]) * TS.dti - GMV.V.tendencies[k] += (GMV.V.new[k] - GMV.V.values[k]) * TS.dti +function GMV_update(grid::Grid, GMV::GridMeanVariables, TS::TimeStepping) + @inbounds for k in real_center_indices(grid) + GMV.U.values[k] += GMV.U.tendencies[k] * TS.dt + GMV.V.values[k] += GMV.V.tendencies[k] * TS.dt + GMV.H.values[k] += GMV.H.tendencies[k] * TS.dt + GMV.QT.values[k] += GMV.QT.tendencies[k] * TS.dt end - return end + function compute_eddy_diffusivities_tke(edmf::EDMF_PrognosticTKE, GMV::GridMeanVariables, Case::CasesBase) grid = get_grid(edmf) param_set = parameter_set(GMV) diff --git a/src/Variables.jl b/src/Variables.jl index 829b5610eb..4d0b1ca014 100644 --- a/src/Variables.jl +++ b/src/Variables.jl @@ -1,14 +1,3 @@ -function update(self::GridMeanVariables, TS::TimeStepping) - grid = self.grid - @inbounds for k in real_center_indices(grid) - self.U.values[k] += self.U.tendencies[k] * TS.dt - self.V.values[k] += self.V.tendencies[k] * TS.dt - self.H.values[k] += self.H.tendencies[k] * TS.dt - self.QT.values[k] += self.QT.tendencies[k] * TS.dt - end - return -end - function initialize_io(self::GridMeanVariables, Stats::NetCDFIO_Stats) add_profile(Stats, "u_mean") add_profile(Stats, "v_mean") diff --git a/src/types.jl b/src/types.jl index 890fe41719..947cf089d6 100644 --- a/src/types.jl +++ b/src/types.jl @@ -227,7 +227,6 @@ end struct VariablePrognostic{T} values::T - new::T tendencies::T loc::String bc::String @@ -238,12 +237,11 @@ struct VariablePrognostic{T} # Value at the current timestep values = field(grid, loc) # Value at the next timestep, used for calculating turbulence tendencies - new = field(grid, loc) tendencies = field(grid, loc) if kind != "scalar" && kind != "velocity" print("Invalid kind setting for variable! Must be scalar or velocity") end - return new{typeof(values)}(values, new, tendencies, loc, bc, kind, name, units) + return new{typeof(values)}(values, tendencies, loc, bc, kind, name, units) end end @@ -1021,10 +1019,10 @@ mutable struct EDMF_PrognosticTKE{PS, A1, A2, IE} # Vertical fluxes for output massflux_h = face_field(grid) massflux_qt = face_field(grid) - diffusive_flux_h = center_field(grid) - diffusive_flux_qt = center_field(grid) - diffusive_flux_u = center_field(grid) - diffusive_flux_v = center_field(grid) + diffusive_flux_h = face_field(grid) + diffusive_flux_qt = face_field(grid) + diffusive_flux_u = face_field(grid) + diffusive_flux_v = face_field(grid) massflux_tke = center_field(grid) # Initialize SDE parameters