explict gmv equations

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"
 

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())

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)

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")

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