Merge pull request #3653 from CliMA/cc/implicit_up_mf_flag

Add flag for running updraft massflux implicitly in GM equations

Author: costachris

config/default_configs/default_config.yml

@@ -315,6 +315,9 @@ zero_tendency:
 implicit_sgs_advection:
   help: "Whether to treat the subgrid-scale vertical advection tendency implicitly [`false` (default), `true`]"
   value: false
+implicit_sgs_mass_flux:
+  help: "Whether to treat the subgrid-scale mass flux tendency implicitly or explicitly in grid-mean equations. Currently updraft only with Jacobian terms 0. [`false` (default), `true`]"
+  value: false
 edmf_coriolis:
   help: "EDMF coriolis [`nothing` (default), `Bomex`,`LifeCycleTan2018`,`Rico`,`ARM_SGP`,`DYCOMS_RF01`,`DYCOMS_RF02`,`GABLS`]"
   value: ~

config/model_configs/prognostic_edmfx_soares_column.yml

@@ -3,7 +3,8 @@ surface_setup: "Soares"
 rayleigh_sponge: true
 turbconv: "prognostic_edmfx"
 implicit_diffusion: true
-implicit_sgs_advection: false
+implicit_sgs_advection: true
+implicit_sgs_mass_flux: true
 approximate_linear_solve_iters: 2
 edmfx_upwinding: "first_order"
 edmfx_entr_model: "Generalized"

src/prognostic_equations/implicit/implicit_solver.jl

@@ -98,6 +98,7 @@ struct ImplicitEquationJacobian{
     F1 <: DerivativeFlag,
     F2 <: DerivativeFlag,
     F3 <: DerivativeFlag,
+    F4 <: DerivativeFlag,
     T <: Fields.FieldVector,
     R <: Base.RefValue,
 }
@@ -111,6 +112,7 @@ struct ImplicitEquationJacobian{
     diffusion_flag::F1
     topography_flag::F2
     sgs_advection_flag::F3
+    sgs_mass_flux_flag::F4
 
     # required by Krylov.jl to evaluate ldiv! with AbstractVector inputs
     temp_b::T
@@ -128,6 +130,7 @@ function Base.zero(jac::ImplicitEquationJacobian)
         jac.diffusion_flag,
         jac.topography_flag,
         jac.sgs_advection_flag,
+        jac.sgs_mass_flux_flag,
         jac.temp_b,
         jac.temp_x,
         jac.transform_flag,
@@ -145,6 +148,8 @@ function ImplicitEquationJacobian(
                       IgnoreDerivative(),
     sgs_advection_flag = atmos.sgs_adv_mode == Implicit() ? UseDerivative() :
                          IgnoreDerivative(),
+    sgs_mass_flux_flag = atmos.sgs_mf_mode == Implicit() ? UseDerivative() :
+                         IgnoreDerivative(),
     transform_flag = false,
 )
     FT = Spaces.undertype(axes(Y.c))
@@ -297,11 +302,28 @@ function ImplicitEquationJacobian(
         ()
     end
 
+    sgs_massflux_blocks = if atmos.turbconv_model isa PrognosticEDMFX
+        @assert n_prognostic_mass_flux_subdomains(atmos.turbconv_model) == 1
+        if use_derivative(sgs_mass_flux_flag)
+            (
+                (@name(c.ρe_tot), @name(c.sgsʲs.:(1).mse)) =>
+                    similar(Y.c, TridiagonalRow),
+                (@name(c.ρq_tot), @name(c.sgsʲs.:(1).q_tot)) =>
+                    similar(Y.c, TridiagonalRow),
+            )
+        else
+            ()
+        end
+    else
+        ()
+    end
+
     matrix = MatrixFields.FieldMatrix(
         identity_blocks...,
         sgs_advection_blocks...,
         advection_blocks...,
         diffusion_blocks...,
+        sgs_massflux_blocks...,
     )
 
     sgs_names_if_available = if atmos.turbconv_model isa PrognosticEDMFX
@@ -381,6 +403,7 @@ function ImplicitEquationJacobian(
         diffusion_flag,
         topography_flag,
         sgs_advection_flag,
+        sgs_mass_flux_flag,
         similar(Y),
         similar(Y),
         transform_flag,
@@ -498,7 +521,13 @@ end
 
 function update_implicit_equation_jacobian!(A, Y, p, dtγ)
     dtγ = float(dtγ)
-    (; matrix, diffusion_flag, sgs_advection_flag, topography_flag) = A
+    (;
+        matrix,
+        diffusion_flag,
+        sgs_advection_flag,
+        topography_flag,
+        sgs_mass_flux_flag,
+    ) = A
     (; ᶜspecific, ᶜK, ᶜts, ᶜp, ᶜΦ, ᶠgradᵥ_ᶜΦ, ᶜh_tot) = p
     (;
         ᶜtemp_C3,
@@ -945,6 +974,42 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
                     dtγ * ᶠtridiagonal_matrix_c3 ⋅
                     DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃))) - (I_u₃,)
             end
+            # add updraft mass flux contributions to grid-mean
+            if use_derivative(sgs_mass_flux_flag)
+
+                (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = p
+                (; bdmr_l, bdmr_r, bdmr) = p
+                is_third_order = edmfx_upwinding == Val(:third_order)
+                ᶠupwind = is_third_order ? ᶠupwind3 : ᶠupwind1
+                ᶠset_upwind_bcs = Operators.SetBoundaryOperator(;
+                    top = Operators.SetValue(zero(CT3{FT})),
+                    bottom = Operators.SetValue(zero(CT3{FT})),
+                ) # Need to wrap ᶠupwind in this for well-defined boundaries.
+                UpwindMatrixRowType =
+                    is_third_order ? QuaddiagonalMatrixRow : BidiagonalMatrixRow
+                ᶠupwind_matrix =
+                    is_third_order ? ᶠupwind3_matrix : ᶠupwind1_matrix
+
+                ᶜkappa_mʲ = p.ᶜtemp_scalar
+                @. ᶜkappa_mʲ =
+                    TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1)) /
+                    TD.cv_m(thermo_params, ᶜtsʲs.:(1))
+
+                # Placeholders for jacobian contributions of updraft massflux to grid-mean
+                # Initialize all derivatives to zero
+
+                ## grid-mean ρρe_tot
+                ∂ᶜρe_tot_err_∂ᶜmseʲ =
+                    matrix[@name(c.ρe_tot), @name(c.sgsʲs.:(1).mse)]
+                ∂ᶜρe_tot_err_∂ᶜmseʲ .= (zero(eltype(∂ᶜρe_tot_err_∂ᶜmseʲ)),)
+
+                ## grid-mean ρq_tot
+                ∂ᶜρq_tot_err_∂ᶜq_totʲ =
+                    matrix[@name(c.ρq_tot), @name(c.sgsʲs.:(1).q_tot)]
+                ∂ᶜρq_tot_err_∂ᶜq_totʲ .= (zero(eltype(∂ᶜρq_tot_err_∂ᶜq_totʲ)),)
+
+            end
+
         elseif rs isa RayleighSponge
             ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ =
                 matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)]
@@ -955,5 +1020,4 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ)
                 ) - (I_u₃,)
         end
     end
-
 end

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -29,6 +29,11 @@ NVTX.@annotate function implicit_tendency!(Yₜ, Y, p, t)
         )
         edmfx_sgs_diffusive_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     end
+
+    if p.atmos.sgs_mf_mode == Implicit()
+        edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    end
+
     # NOTE: All ρa tendencies should be applied before calling this function
     pressure_work_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
 

src/prognostic_equations/remaining_tendency.jl

@@ -111,7 +111,9 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
 
     radiation_tendency!(Yₜ, Y, p, t, p.atmos.radiation_mode)
     edmfx_entr_detr_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
-    edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    if p.atmos.sgs_mf_mode == Explicit()
+        edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
+    end
     edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     edmfx_filter_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)
     edmfx_tke_tendency!(Yₜ, Y, p, t, p.atmos.turbconv_model)

src/solver/type_getters.jl

@@ -55,6 +55,9 @@ function get_atmos(config::AtmosConfig, params)
     implicit_sgs_advection = parsed_args["implicit_sgs_advection"]
     @assert implicit_sgs_advection in (true, false)
 
+    implicit_sgs_mass_flux = parsed_args["implicit_sgs_mass_flux"]
+    @assert implicit_sgs_mass_flux in (true, false)
+
     edmfx_model = EDMFXModel(;
         entr_model = get_entrainment_model(parsed_args),
         detr_model = get_detrainment_model(parsed_args),
@@ -96,6 +99,7 @@ function get_atmos(config::AtmosConfig, params)
         vert_diff,
         diff_mode = implicit_diffusion ? Implicit() : Explicit(),
         sgs_adv_mode = implicit_sgs_advection ? Implicit() : Explicit(),
+        sgs_mf_mode = implicit_sgs_mass_flux ? Implicit() : Explicit(),
         viscous_sponge = get_viscous_sponge_model(parsed_args, params, FT),
         smagorinsky_lilly = get_smagorinsky_lilly_model(parsed_args),
         rayleigh_sponge = get_rayleigh_sponge_model(parsed_args, params, FT),

src/solver/types.jl

@@ -516,6 +516,7 @@ Base.@kwdef struct AtmosModel{
     VD,
     DM,
     SAM,
+    SMM,
     VS,
     SL,
     RS,
@@ -553,6 +554,7 @@ Base.@kwdef struct AtmosModel{
     vert_diff::VD = nothing
     diff_mode::DM = nothing
     sgs_adv_mode::SAM = nothing
+    sgs_mf_mode::SMM = nothing
     viscous_sponge::VS = nothing
     smagorinsky_lilly::SL = nothing
     rayleigh_sponge::RS = nothing