From 6c29cf059c1569573081bccc98b022edebdea97c Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Mon, 29 Jan 2024 17:08:16 -0800 Subject: [PATCH] make updraft advection and diffusion implicit in prognostic edmf --- .../gpu_prognostic_edmfx_aquaplanet.yml | 2 +- .../prognostic_edmfx_aquaplanet.yml | 2 +- .../prognostic_edmfx_bomex_box.yml | 6 +- .../prognostic_edmfx_dycoms_rf01_box.yml | 6 +- .../prognostic_edmfx_gabls_box.yml | 8 +- .../prognostic_edmfx_simpleplume_column.yml | 11 +- .../prognostic_edmfx_trmm_column.yml | 3 +- ...erf_target_prognostic_edmfx_aquaplanet.yml | 5 +- src/cache/cache.jl | 2 + src/prognostic_equations/advection.jl | 6 +- .../implicit/implicit_solver.jl | 297 +++++++++++++++++- .../implicit/implicit_tendency.jl | 8 + .../remaining_tendency.jl | 16 +- src/surface_conditions/surface_setups.jl | 2 +- src/utils/abbreviations.jl | 4 + toml/prognostic_edmfx_bomex_box.toml | 2 +- toml/prognostic_edmfx_box.toml | 12 +- toml/prognostic_edmfx_dycoms_rf01_box.toml | 12 +- 18 files changed, 361 insertions(+), 43 deletions(-) diff --git a/config/gpu_configs/gpu_prognostic_edmfx_aquaplanet.yml b/config/gpu_configs/gpu_prognostic_edmfx_aquaplanet.yml index fee1ea6dc2..43972f8099 100644 --- a/config/gpu_configs/gpu_prognostic_edmfx_aquaplanet.yml +++ b/config/gpu_configs/gpu_prognostic_edmfx_aquaplanet.yml @@ -5,7 +5,7 @@ turbconv: prognostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "Generalized" -edmfx_detr_model: "ConstantArea" +edmfx_detr_model: "Generalized" edmfx_nh_pressure: true edmfx_velocity_relaxation: true edmfx_sgs_mass_flux: true diff --git a/config/model_configs/prognostic_edmfx_aquaplanet.yml b/config/model_configs/prognostic_edmfx_aquaplanet.yml index 177455a8d9..cb429c0a73 100644 --- a/config/model_configs/prognostic_edmfx_aquaplanet.yml +++ b/config/model_configs/prognostic_edmfx_aquaplanet.yml @@ -5,7 +5,7 @@ turbconv: prognostic_edmfx prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "Generalized" -edmfx_detr_model: "ConstantArea" +edmfx_detr_model: "Generalized" edmfx_nh_pressure: true edmfx_velocity_relaxation: true edmfx_sgs_mass_flux: true diff --git a/config/model_configs/prognostic_edmfx_bomex_box.yml b/config/model_configs/prognostic_edmfx_bomex_box.yml index 8d670d5c21..6a9e419895 100644 --- a/config/model_configs/prognostic_edmfx_bomex_box.yml +++ b/config/model_configs/prognostic_edmfx_bomex_box.yml @@ -5,12 +5,16 @@ edmf_coriolis: "Bomex" ls_adv: "Bomex" surface_setup: "Bomex" turbconv: "prognostic_edmfx" +implicit_diffusion: true +approximate_linear_solve_iters: 2 +max_newton_iters_ode: 3 edmfx_upwinding: first_order edmfx_entr_model: "Generalized" edmfx_detr_model: "Generalized" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true +edmfx_velocity_relaxation: true prognostic_tke: true moist: "equil" config: "box" @@ -23,7 +27,7 @@ y_elem: 2 z_elem: 60 z_stretch: false perturb_initstate: false -dt: "10secs" +dt: "50secs" t_end: "6hours" dt_save_state_to_disk: "30mins" toml: [toml/prognostic_edmfx_bomex_box.toml] diff --git a/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml b/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml index f6a68736a4..57cdd14155 100644 --- a/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml +++ b/config/model_configs/prognostic_edmfx_dycoms_rf01_box.yml @@ -5,6 +5,9 @@ edmf_coriolis: DYCOMS_RF01 rad: DYCOMS_RF01 surface_setup: DYCOMS_RF01 turbconv: prognostic_edmfx +implicit_diffusion: true +approximate_linear_solve_iters: 2 +max_newton_iters_ode: 3 edmfx_upwinding: first_order edmfx_sgsflux_upwinding: first_order edmfx_entr_model: "Generalized" @@ -12,6 +15,7 @@ edmfx_detr_model: "Generalized" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true +edmfx_velocity_relaxation: true prognostic_tke: true moist: equil config: box @@ -23,7 +27,7 @@ x_elem: 2 y_elem: 2 z_elem: 30 z_stretch: false -dt: 10secs +dt: 40secs t_end: 4hours dt_save_state_to_disk: 10mins toml: [toml/prognostic_edmfx_dycoms_rf01_box.toml] diff --git a/config/model_configs/prognostic_edmfx_gabls_box.yml b/config/model_configs/prognostic_edmfx_gabls_box.yml index 943a676a1c..db7891b0bc 100644 --- a/config/model_configs/prognostic_edmfx_gabls_box.yml +++ b/config/model_configs/prognostic_edmfx_gabls_box.yml @@ -3,12 +3,16 @@ initial_condition: GABLS edmf_coriolis: GABLS surface_setup: GABLS turbconv: "prognostic_edmfx" +implicit_diffusion: true +approximate_linear_solve_iters: 2 +max_newton_iters_ode: 3 edmfx_upwinding: "first_order" edmfx_entr_model: "Generalized" -edmfx_detr_model: "ConstantArea" +edmfx_detr_model: "Generalized" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true +edmfx_velocity_relaxation: true prognostic_tke: true moist: "equil" config: "box" @@ -20,7 +24,7 @@ x_elem: 2 y_elem: 2 z_elem: 8 z_stretch: false -dt: "100secs" +dt: "60secs" t_end: "9hours" dt_save_state_to_disk: "30mins" perturb_initstate: false diff --git a/config/model_configs/prognostic_edmfx_simpleplume_column.yml b/config/model_configs/prognostic_edmfx_simpleplume_column.yml index 9625230232..6d387684d0 100644 --- a/config/model_configs/prognostic_edmfx_simpleplume_column.yml +++ b/config/model_configs/prognostic_edmfx_simpleplume_column.yml @@ -2,22 +2,25 @@ job_id: "prognostic_edmfx_simpleplume_column" initial_condition: "SimplePlume" surface_setup: "SimplePlume" turbconv: "prognostic_edmfx" +implicit_diffusion: true +approximate_linear_solve_iters: 2 +max_newton_iters_ode: 3 gs_tendency: false edmfx_upwinding: first_order edmfx_entr_model: "Generalized" edmfx_detr_model: "Generalized" -edmfx_sgs_mass_flux: true -edmfx_sgs_diffusive_flux: true +edmfx_sgs_mass_flux: false +edmfx_sgs_diffusive_flux: false edmfx_nh_pressure: true edmfx_velocity_relaxation: true -prognostic_tke: true +prognostic_tke: false moist: "equil" config: "column" z_max: 4e3 z_elem: 80 z_stretch: false perturb_initstate: false -dt: "1secs" +dt: 20secs" t_end: "12hours" dt_save_to_disk: "10mins" toml: [toml/prognostic_edmfx_simpleplume.toml] diff --git a/config/model_configs/prognostic_edmfx_trmm_column.yml b/config/model_configs/prognostic_edmfx_trmm_column.yml index 8b67b181a8..7bb47d44fb 100644 --- a/config/model_configs/prognostic_edmfx_trmm_column.yml +++ b/config/model_configs/prognostic_edmfx_trmm_column.yml @@ -5,10 +5,11 @@ surface_setup: TRMM_LBA turbconv: prognostic_edmfx edmfx_upwinding: first_order edmfx_entr_model: "Generalized" -edmfx_detr_model: "ConstantArea" +edmfx_detr_model: "Generalized" edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true edmfx_nh_pressure: true +edmfx_velocity_relaxation: true prognostic_tke: true moist: equil apply_limiter: false diff --git a/config/perf_configs/flame_perf_target_prognostic_edmfx_aquaplanet.yml b/config/perf_configs/flame_perf_target_prognostic_edmfx_aquaplanet.yml index 06cfebdc5a..ab256008bb 100644 --- a/config/perf_configs/flame_perf_target_prognostic_edmfx_aquaplanet.yml +++ b/config/perf_configs/flame_perf_target_prognostic_edmfx_aquaplanet.yml @@ -1,12 +1,13 @@ job_id: flame_perf_target_prognostic_edmfx_aquaplanet -surface_setup: DefaultExchangeCoefficients +surface_setup: DefaultExchangeCoefficients rad: gray vert_diff: "false" turbconv: prognostic_edmfx +implicit_diffusion: false prognostic_tke: true edmfx_upwinding: first_order edmfx_entr_model: "Generalized" -edmfx_detr_model: "ConstantArea" +edmfx_detr_model: "Generalized" edmfx_nh_pressure: true edmfx_sgs_mass_flux: true edmfx_sgs_diffusive_flux: true diff --git a/src/cache/cache.jl b/src/cache/cache.jl index 81d37a5151..1131dde27d 100644 --- a/src/cache/cache.jl +++ b/src/cache/cache.jl @@ -117,6 +117,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info) ᶠcoord = Fields.local_geometry_field(Y.f).coordinates grav = FT(CAP.grav(params)) ᶜΦ = grav .* ᶜcoord.z + ᶠΦ = grav .* ᶠcoord.z if atmos.numerics.use_reference_state R_d = FT(CAP.R_d(params)) @@ -186,6 +187,7 @@ function build_cache(Y, atmos, params, surface_setup, sim_info) core = ( ᶜΦ, ᶠgradᵥ_ᶜΦ = ᶠgradᵥ.(ᶜΦ), + ᶜgradᵥ_ᶠΦ = ᶜgradᵥ.(ᶠΦ), ᶜρ_ref, ᶜp_ref, ᶜT = similar(Y.c, FT), diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 06c0d44bfa..5fc5f16dba 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -217,15 +217,13 @@ function edmfx_sgs_vertical_advection_tendency!( ᶜa_scalar = p.scratch.ᶜtemp_scalar ᶜu₃ʲ = p.scratch.ᶜtemp_C3 ᶜKᵥʲ = p.scratch.ᶜtemp_scalar_2 - ᶜinterp_lb = Operators.LeftBiasedF2C() - ᶜinterp_rb = Operators.RightBiasedF2C() for j in 1:n # TODO: Add a biased GradientF2F operator in ClimaCore @. ᶜu₃ʲ[colidx] = ᶜinterp(Y.f.sgsʲs.:($$j).u₃[colidx]) @. ᶜKᵥʲ[colidx] = ifelse( ᶜu₃ʲ[colidx].components.data.:1 > 0, - ᶜinterp_lb(ᶠKᵥʲs.:($$j)[colidx]), - ᶜinterp_rb(ᶠKᵥʲs.:($$j)[colidx]), + ᶜleft_bias(ᶠKᵥʲs.:($$j)[colidx]), + ᶜright_bias(ᶠKᵥʲs.:($$j)[colidx]), ) # For the updraft u_3 equation, we assume the grid-mean to be hydrostatic # and calcuate the buoyancy term relative to the grid-mean density. diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl index 0b126ddcef..3c2af64020 100644 --- a/src/prognostic_equations/implicit/implicit_solver.jl +++ b/src/prognostic_equations/implicit/implicit_solver.jl @@ -125,12 +125,15 @@ function ImplicitEquationJacobian( FT = Spaces.undertype(axes(Y.c)) CTh = CTh_vector_type(axes(Y.c)) + DiagonalRow = DiagonalMatrixRow{FT} TridiagonalRow = TridiagonalMatrixRow{FT} BidiagonalRow_C3 = BidiagonalMatrixRow{C3{FT}} TridiagonalRow_ACTh = TridiagonalMatrixRow{Adjoint{FT, CTh{FT}}} BidiagonalRow_ACT3 = BidiagonalMatrixRow{Adjoint{FT, CT3{FT}}} BidiagonalRow_C3xACTh = BidiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CTh{FT})')} + DiagonalRow_C3xACT3 = + DiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')} TridiagonalRow_C3xACT3 = TridiagonalMatrixRow{typeof(zero(C3{FT}) * zero(CT3{FT})')} @@ -139,6 +142,7 @@ function ImplicitEquationJacobian( ρq_tot_if_available = is_in_Y(@name(c.ρq_tot)) ? (@name(c.ρq_tot),) : () ρatke_if_available = is_in_Y(@name(c.sgs⁰.ρatke)) ? (@name(c.sgs⁰.ρatke),) : () + sfc_if_available = is_in_Y(@name(sfc)) ? (@name(sfc),) : () tracer_names = ( @name(c.ρq_tot), @@ -148,14 +152,12 @@ function ImplicitEquationJacobian( @name(c.ρq_sno), ) available_tracer_names = MatrixFields.unrolled_filter(is_in_Y, tracer_names) - other_names = (@name(c.sgsʲs), @name(f.sgsʲs), @name(sfc)) - other_available_names = MatrixFields.unrolled_filter(is_in_Y, other_names) # Note: We have to use FT(-1) * I instead of -I because inv(-1) == -1.0, # which means that multiplying inv(-1) by a Float32 will yield a Float64. identity_blocks = MatrixFields.unrolled_map( name -> (name, name) => FT(-1) * I, - (@name(c.ρ), other_available_names...), + (@name(c.ρ), sfc_if_available...), ) active_scalar_names = (@name(c.ρ), @name(c.ρe_tot), ρq_tot_if_available...) @@ -212,24 +214,102 @@ function ImplicitEquationJacobian( ) end + sgs_blocks = if atmos.turbconv_model isa PrognosticEDMFX + @assert n_prognostic_mass_flux_subdomains(atmos.turbconv_model) == 1 + sgs_scalar_names = ( + @name(c.sgsʲs.:(1).q_tot), + @name(c.sgsʲs.:(1).mse), + @name(c.sgsʲs.:(1).ρa) + ) + if use_derivative(diffusion_flag) + ( + MatrixFields.unrolled_map( + name -> (name, name) => similar(Y.c, TridiagonalRow), + sgs_scalar_names, + )..., + (@name(c.sgsʲs.:(1).mse), @name(c.ρ)) => + similar(Y.c, DiagonalRow), + (@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)) => + similar(Y.c, DiagonalRow), + (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)) => + similar(Y.c, TridiagonalRow), + (@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)) => + similar(Y.c, TridiagonalRow), + (@name(f.sgsʲs.:(1).u₃), @name(c.ρ)) => + similar(Y.f, BidiagonalRow_C3), + (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)) => + similar(Y.f, BidiagonalRow_C3), + (@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)) => + similar(Y.f, BidiagonalRow_C3), + (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => + similar(Y.f, TridiagonalRow_C3xACT3), + ) + else + ( + MatrixFields.unrolled_map( + name -> (name, name) => FT(-1) * I, + sgs_scalar_names, + )..., + (@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)) => + !isnothing(atmos.rayleigh_sponge) ? + similar(Y.f, DiagonalRow_C3xACT3) : FT(-1) * I, + ) + end + else + () + end + matrix = MatrixFields.FieldMatrix( identity_blocks..., + sgs_blocks..., advection_blocks..., diffusion_blocks..., ) - names₁_group₁ = (@name(c.ρ), other_available_names...) + sgs_names_if_available = if atmos.turbconv_model isa PrognosticEDMFX + ( + @name(c.sgsʲs.:(1).q_tot), + @name(c.sgsʲs.:(1).mse), + @name(c.sgsʲs.:(1).ρa), + @name(f.sgsʲs.:(1).u₃), + ) + else + () + end + names₁_group₁ = (@name(c.ρ), sfc_if_available...) names₁_group₂ = (available_tracer_names..., ρatke_if_available...) names₁_group₃ = (@name(c.ρe_tot),) - names₁ = (names₁_group₁..., names₁_group₂..., names₁_group₃...) + names₁ = ( + names₁_group₁..., + sgs_names_if_available..., + names₁_group₂..., + names₁_group₃..., + ) alg₂ = MatrixFields.BlockLowerTriangularSolve(@name(c.uₕ)) alg = if use_derivative(diffusion_flag) - alg₁_subalg₂ = + alg₁_subalg₂ = if atmos.turbconv_model isa PrognosticEDMFX + (; + alg₂ = MatrixFields.BlockLowerTriangularSolve( + @name(c.sgsʲs.:(1).q_tot); + alg₂ = MatrixFields.BlockLowerTriangularSolve( + @name(c.sgsʲs.:(1).mse); + alg₂ = MatrixFields.BlockLowerTriangularSolve( + @name(c.sgsʲs.:(1).ρa), + @name(f.sgsʲs.:(1).u₃); + alg₂ = MatrixFields.BlockLowerTriangularSolve( + names₁_group₂..., + ), + ), + ), + ) + ) + else is_in_Y(@name(c.ρq_tot)) ? (; alg₂ = MatrixFields.BlockLowerTriangularSolve(names₁_group₂...) ) : (;) + end alg₁ = MatrixFields.BlockLowerTriangularSolve( names₁_group₁...; alg₁_subalg₂..., @@ -317,13 +397,20 @@ NVTX.@annotate function Wfact!(A, Y, p, dtγ, t) ( use_derivative(A.diffusion_flag) && p.atmos.turbconv_model isa PrognosticEDMFX ? - (; p.precomputed.ᶜρa⁰) : (;) + (; + p.core.ᶜgradᵥ_ᶠΦ, + p.precomputed.ᶜρa⁰, + p.precomputed.ᶜρʲs, + p.precomputed.ᶠu³ʲs, + p.precomputed.ᶜtsʲs, + ) : (;) )..., p.core.ᶜΦ, p.core.ᶠgradᵥ_ᶜΦ, p.core.ᶜρ_ref, p.core.ᶜp_ref, p.scratch.ᶜtemp_scalar, + p.scratch.ᶜtemp_C3, p.scratch.ᶠtemp_CT3, p.scratch.∂ᶜK_∂ᶜuₕ, p.scratch.∂ᶜK_∂ᶠu₃, @@ -355,6 +442,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) (; ᶜspecific, ᶜK, ᶜts, ᶜp, ᶜΦ, ᶠgradᵥ_ᶜΦ, ᶜρ_ref, ᶜp_ref) = p (; ∂ᶜK_∂ᶜuₕ, ∂ᶜK_∂ᶠu₃, ᶠp_grad_matrix, ᶜadvection_matrix) = p (; ᶜdiffusion_h_matrix, ᶜdiffusion_u_matrix, params) = p + (; edmfx_upwinding) = p.atmos.numerics FT = Spaces.undertype(axes(Y.c)) CTh = CTh_vector_type(axes(Y.c)) @@ -593,4 +681,199 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) DiagonalMatrixRow(-(ᶜwₚ[colidx]) / ᶜρ[colidx]) end end + + if p.atmos.turbconv_model isa PrognosticEDMFX + if use_derivative(diffusion_flag) + (; ᶜgradᵥ_ᶠΦ, ᶜρʲs, ᶠu³ʲs, ᶜtsʲs) = 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 + ᶠset_upwind_matrix_bcs = Operators.SetBoundaryOperator(; + top = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})), + bottom = Operators.SetValue(zero(UpwindMatrixRowType{CT3{FT}})), + ) # Need to wrap ᶠupwind_matrix in this for well-defined boundaries. + ᶜkappa_mʲ = p.ᶜtemp_scalar + @. ᶜkappa_mʲ[colidx] = + TD.gas_constant_air(thermo_params, ᶜtsʲs.:(1)[colidx]) / + TD.cv_m(thermo_params, ᶜtsʲs.:(1)[colidx]) + ∂ᶜq_totʲ_err_∂ᶜq_totʲ = + matrix[@name(c.sgsʲs.:(1).q_tot), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶜq_totʲ_err_∂ᶜq_totʲ[colidx] = + dtγ * ( + DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1)[colidx])) - + ᶜadvdivᵥ_matrix() ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) + ) - (I,) + + ∂ᶜmseʲ_err_∂ᶜq_totʲ = + matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶜmseʲ_err_∂ᶜq_totʲ[colidx] = + dtγ * ( + -DiagonalMatrixRow( + adjoint(ᶜinterp(ᶠu³ʲs.:(1)[colidx])) * + ᶜgradᵥ_ᶠΦ[colidx] * + Y.c.ρ[colidx] * + ᶜkappa_mʲ[colidx] / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0, + ) + ) + ∂ᶜmseʲ_err_∂ᶜρ = matrix[@name(c.sgsʲs.:(1).mse), @name(c.ρ)] + @. ∂ᶜmseʲ_err_∂ᶜρ[colidx] = + dtγ * ( + -DiagonalMatrixRow( + adjoint(ᶜinterp(ᶠu³ʲs.:(1)[colidx])) * + ᶜgradᵥ_ᶠΦ[colidx] / ᶜρʲs.:(1)[colidx], + ) + ) + ∂ᶜmseʲ_err_∂ᶜmseʲ = + matrix[@name(c.sgsʲs.:(1).mse), @name(c.sgsʲs.:(1).mse)] + @. ∂ᶜmseʲ_err_∂ᶜmseʲ[colidx] = + dtγ * ( + DiagonalMatrixRow(ᶜadvdivᵥ(ᶠu³ʲs.:(1)[colidx])) - + ᶜadvdivᵥ_matrix() ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) - + DiagonalMatrixRow( + adjoint(ᶜinterp(ᶠu³ʲs.:(1)[colidx])) * + ᶜgradᵥ_ᶠΦ[colidx] * + Y.c.ρ[colidx] * + ᶜkappa_mʲ[colidx] / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]), + ) + ) - (I,) + + ∂ᶜρaʲ_err_∂ᶜq_totʲ = + matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶜρaʲ_err_∂ᶜq_totʲ[colidx] = + dtγ * ᶜadvdivᵥ_matrix() ⋅ ( + DiagonalMatrixRow( + ᶠset_upwind_bcs( + ᶠupwind( + ᶠu³ʲs.:(1)[colidx], + draft_area( + Y.c.sgsʲs.:(1).ρa[colidx], + ᶜρʲs.:(1)[colidx], + ), + ), + ), + ) ⋅ ᶠwinterp_matrix(ᶜJ[colidx]) ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0, + ) - + DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅ + DiagonalMatrixRow( + Y.c.sgsʲs.:(1).ρa[colidx] * ᶜkappa_mʲ[colidx] / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0, + ) + ) + ∂ᶜρaʲ_err_∂ᶜmseʲ = + matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).mse)] + @. ∂ᶜρaʲ_err_∂ᶜmseʲ[colidx] = + dtγ * ᶜadvdivᵥ_matrix() ⋅ ( + DiagonalMatrixRow( + ᶠset_upwind_bcs( + ᶠupwind( + ᶠu³ʲs.:(1)[colidx], + draft_area( + Y.c.sgsʲs.:(1).ρa[colidx], + ᶜρʲs.:(1)[colidx], + ), + ), + ), + ) ⋅ ᶠwinterp_matrix(ᶜJ[colidx]) ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]), + ) - + DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅ + DiagonalMatrixRow( + Y.c.sgsʲs.:(1).ρa[colidx] * ᶜkappa_mʲ[colidx] / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]), + ) + ) + ∂ᶜρaʲ_err_∂ᶜρaʲ = + matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)] + @. ∂ᶜρaʲ_err_∂ᶜρaʲ[colidx] = + dtγ * ᶜadvdivᵥ_matrix() ⋅ ( + -DiagonalMatrixRow( + ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx]), + ) ⋅ + ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅ + DiagonalMatrixRow(1 / ᶜρʲs.:(1)[colidx]) + ) - (I,) + + ∂ᶠu₃ʲ_err_∂ᶜρ = matrix[@name(f.sgsʲs.:(1).u₃), @name(c.ρ)] + @. ∂ᶠu₃ʲ_err_∂ᶜρ[colidx] = + dtγ * DiagonalMatrixRow( + ᶠgradᵥ_ᶜΦ[colidx] / ᶠinterp(ᶜρʲs.:(1)[colidx]), + ) ⋅ ᶠinterp_matrix() + ∂ᶠu₃ʲ_err_∂ᶜq_totʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).q_tot)] + @. ∂ᶠu₃ʲ_err_∂ᶜq_totʲ[colidx] = + dtγ * DiagonalMatrixRow( + ᶠgradᵥ_ᶜΦ[colidx] * ᶠinterp(Y.c.ρ[colidx]) / + (ᶠinterp(ᶜρʲs.:(1)[colidx]))^2, + ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0, + ) + ∂ᶠu₃ʲ_err_∂ᶜmseʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(c.sgsʲs.:(1).mse)] + @. ∂ᶠu₃ʲ_err_∂ᶜmseʲ[colidx] = + dtγ * DiagonalMatrixRow( + ᶠgradᵥ_ᶜΦ[colidx] * ᶠinterp(Y.c.ρ[colidx]) / + (ᶠinterp(ᶜρʲs.:(1)[colidx]))^2, + ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( + ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / + ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]), + ) + ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)] + ᶜu₃ʲ = p.ᶜtemp_C3 + @. ᶜu₃ʲ[colidx] = ᶜinterp(Y.f.sgsʲs.:(1).u₃[colidx]) + if p.atmos.rayleigh_sponge isa RayleighSponge + @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ[colidx] = + dtγ * -( + ᶠgradᵥ_matrix() ⋅ ifelse( + ᶜu₃ʲ[colidx].components.data.:1 > 0, + convert( + BidiagonalMatrixRow{FT}, + ᶜleft_bias_matrix(), + ), + convert( + BidiagonalMatrixRow{FT}, + ᶜright_bias_matrix(), + ), + ) ⋅ DiagonalMatrixRow( + adjoint(CT3(Y.f.sgsʲs.:(1).u₃[colidx])), + ) + DiagonalMatrixRow( + p.ᶠβ_rayleigh_w[colidx] * (one_C3xACT3,), + ) + ) - (I_u₃,) + else + @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ[colidx] = + dtγ * -(ᶠgradᵥ_matrix()) ⋅ ifelse( + ᶜu₃ʲ[colidx].components.data.:1 > 0, + convert(BidiagonalMatrixRow{FT}, ᶜleft_bias_matrix()), + convert(BidiagonalMatrixRow{FT}, ᶜright_bias_matrix()), + ) ⋅ + DiagonalMatrixRow(adjoint(CT3(Y.f.sgsʲs.:(1).u₃[colidx]))) - + (I_u₃,) + end + elseif p.atmos.rayleigh_sponge isa RayleighSponge + ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ = + matrix[@name(f.sgsʲs.:(1).u₃), @name(f.sgsʲs.:(1).u₃)] + @. ∂ᶠu₃ʲ_err_∂ᶠu₃ʲ[colidx] = + dtγ * + -DiagonalMatrixRow(p.ᶠβ_rayleigh_w[colidx] * (one_C3xACT3,)) - + (I_u₃,) + end + end + end diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 51c6382ea9..aedcc00b07 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -10,6 +10,14 @@ NVTX.@annotate function implicit_tendency!(Yₜ, Y, p, t) Fields.bycolumn(axes(Y.c)) do colidx implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) if p.atmos.diff_mode == Implicit() + edmfx_sgs_vertical_advection_tendency!( + Yₜ, + Y, + p, + t, + colidx, + p.atmos.turbconv_model, + ) vertical_diffusion_boundary_layer_tendency!( Yₜ, Y, diff --git a/src/prognostic_equations/remaining_tendency.jl b/src/prognostic_equations/remaining_tendency.jl index 07ba4cba15..dd87abe296 100644 --- a/src/prognostic_equations/remaining_tendency.jl +++ b/src/prognostic_equations/remaining_tendency.jl @@ -35,6 +35,14 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) large_scale_advection_tendency!(Yₜ, Y, p, t, colidx, p.atmos.ls_adv) if p.atmos.diff_mode == Explicit() + edmfx_sgs_vertical_advection_tendency!( + Yₜ, + Y, + p, + t, + colidx, + p.atmos.turbconv_model, + ) vertical_diffusion_boundary_layer_tendency!( Yₜ, Y, @@ -54,14 +62,6 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t) end radiation_tendency!(Yₜ, Y, p, t, colidx, p.atmos.radiation_mode) - edmfx_sgs_vertical_advection_tendency!( - Yₜ, - Y, - p, - t, - colidx, - p.atmos.turbconv_model, - ) edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, p.atmos.turbconv_model) edmfx_sgs_mass_flux_tendency!( Yₜ, diff --git a/src/surface_conditions/surface_setups.jl b/src/surface_conditions/surface_setups.jl index 4e96329b3a..881c805bbf 100644 --- a/src/surface_conditions/surface_setups.jl +++ b/src/surface_conditions/surface_setups.jl @@ -237,7 +237,7 @@ function (::SimplePlume)(params) T = FT(310) p = FT(101500) q_vap = FT(0.02245) - θ_flux = FT(8e-3) + θ_flux = FT(8) q_flux = FT(0) z0 = FT(1e-4) ustar = FT(0.28) diff --git a/src/utils/abbreviations.jl b/src/utils/abbreviations.jl index b15bd0444c..84c3428da8 100644 --- a/src/utils/abbreviations.jl +++ b/src/utils/abbreviations.jl @@ -34,6 +34,8 @@ const ᶜadvdivᵥ = Operators.DivergenceF2C( const ᶜprecipdivᵥ = Operators.DivergenceF2C(top = Operators.SetValue(CT3(0))) const ᶠright_bias = Operators.RightBiasedC2F() # for free outflow in ᶜprecipdivᵥ +const ᶜleft_bias = Operators.LeftBiasedF2C() +const ᶜright_bias = Operators.RightBiasedF2C() # TODO: Implement proper extrapolation instead of simply reusing the first # interior value at the surface. @@ -78,6 +80,8 @@ const ᶠfct_zalesak = Operators.FCTZalesak( ) const ᶜinterp_matrix = MatrixFields.operator_matrix(ᶜinterp) +const ᶜleft_bias_matrix = MatrixFields.operator_matrix(ᶜleft_bias) +const ᶜright_bias_matrix = MatrixFields.operator_matrix(ᶜright_bias) const ᶜdivᵥ_matrix = MatrixFields.operator_matrix(ᶜdivᵥ) const ᶜadvdivᵥ_matrix = MatrixFields.operator_matrix(ᶜadvdivᵥ) const ᶜprecipdivᵥ_matrix = MatrixFields.operator_matrix(ᶜprecipdivᵥ) diff --git a/toml/prognostic_edmfx_bomex_box.toml b/toml/prognostic_edmfx_bomex_box.toml index 7f9370356e..20ffe2d1af 100644 --- a/toml/prognostic_edmfx_bomex_box.toml +++ b/toml/prognostic_edmfx_bomex_box.toml @@ -11,7 +11,7 @@ value = 0.7 value = 0.0002 [entr_coeff] -value = 0.2 +value = 0.3 [min_area_limiter_scale] value = 0 diff --git a/toml/prognostic_edmfx_box.toml b/toml/prognostic_edmfx_box.toml index 9a9cb98ae1..7341d979b5 100644 --- a/toml/prognostic_edmfx_box.toml +++ b/toml/prognostic_edmfx_box.toml @@ -5,13 +5,13 @@ value = 0.1 value = 1.0e-5 [EDMF_max_area] -value = 0.9 +value = 0.7 [entr_inv_tau] -value = 0.002 +value = 0.0001 [entr_coeff] -value = 0 +value = 0.1 [min_area_limiter_scale] value = 0 @@ -25,5 +25,11 @@ value = 0 [detr_buoy_coeff] value = 0 +[detr_vertdiv_coeff] +value = 0 + +[detr_massflux_vertdiv_coeff] +value = 1 + [max_area_limiter_scale] value = 0 diff --git a/toml/prognostic_edmfx_dycoms_rf01_box.toml b/toml/prognostic_edmfx_dycoms_rf01_box.toml index b6b1616bbd..05b6c0cb91 100644 --- a/toml/prognostic_edmfx_dycoms_rf01_box.toml +++ b/toml/prognostic_edmfx_dycoms_rf01_box.toml @@ -8,13 +8,13 @@ value = 1.0e-5 value = 0.7 [entr_inv_tau] -value = 0 +value = 0.0001 [entr_coeff] -value = 0.2 +value = 0.1 [min_area_limiter_scale] -value = 0.01 +value = 0 [min_area_limiter_power] value = 100 @@ -23,16 +23,16 @@ value = 100 value = 0 [detr_coeff] -value = 1e-3 +value = 0 [detr_buoy_coeff] -value = 0.3 +value = 0 [detr_vertdiv_coeff] value = 0 [detr_massflux_vertdiv_coeff] -value = 0.5 +value = 1 [max_area_limiter_scale] value = 0.01