diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 5fc5f16dba..ffa9d84ec1 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -73,6 +73,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) point_type = eltype(Fields.coordinate_field(Y.c)) (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J + ᶠJ = Fields.local_geometry_field(Y.f).J (; ᶜf³, ᶠf¹², ᶜΦ) = p.core (; ᶜu, ᶠu³, ᶜK) = p.precomputed (; edmfx_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing @@ -116,6 +117,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) coeff, Yₜ.c.ρe_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], Y.c.ρ[colidx], ᶠu³[colidx], ᶜh_tot[colidx], @@ -132,6 +134,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) coeff, ᶜρχₜ[colidx], ᶜJ[colidx], + ᶠJ[colidx], Y.c.ρ[colidx], ᶠu³[colidx], ᶜχ[colidx], @@ -184,6 +187,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) vertical_transport!( Yₜ.c.sgs⁰.ρatke[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρ⁰[colidx], ᶠu³⁰[colidx], ᶜa_scalar[colidx], @@ -209,6 +213,7 @@ function edmfx_sgs_vertical_advection_tendency!( n = n_prognostic_mass_flux_subdomains(turbconv_model) (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J + ᶠJ = Fields.local_geometry_field(Y.f).J (; edmfx_upwinding) = p.atmos.numerics (; ᶠu³ʲs, ᶠKᵥʲs, ᶜρʲs) = p.precomputed (; ᶠgradᵥ_ᶜΦ) = p.core @@ -244,6 +249,7 @@ function edmfx_sgs_vertical_advection_tendency!( vertical_transport!( Yₜ.c.sgsʲs.:($j).ρa[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³ʲs.:($j)[colidx], ᶜa_scalar[colidx], diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl index 10dba82f65..f9da02e0d1 100644 --- a/src/prognostic_equations/edmfx_sgs_flux.jl +++ b/src/prognostic_equations/edmfx_sgs_flux.jl @@ -20,6 +20,7 @@ function edmfx_sgs_mass_flux_tendency!( (; ᶜρa⁰, ᶜρ⁰, ᶠu³⁰, ᶜK⁰, ᶜmse⁰, ᶜq_tot⁰) = p.precomputed (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J + ᶠJ = Fields.local_geometry_field(Y.f).J if p.atmos.edmfx_sgs_mass_flux # energy @@ -35,6 +36,7 @@ function edmfx_sgs_mass_flux_tendency!( vertical_transport!( Yₜ.c.ρe_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³_diff_colidx, ᶜa_scalar_colidx, @@ -49,6 +51,7 @@ function edmfx_sgs_mass_flux_tendency!( vertical_transport!( Yₜ.c.ρe_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρ⁰[colidx], ᶠu³_diff_colidx, ᶜa_scalar_colidx, @@ -66,6 +69,7 @@ function edmfx_sgs_mass_flux_tendency!( vertical_transport!( Yₜ.c.ρq_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³_diff_colidx, ᶜa_scalar_colidx, @@ -80,6 +84,7 @@ function edmfx_sgs_mass_flux_tendency!( vertical_transport!( Yₜ.c.ρq_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρ⁰[colidx], ᶠu³_diff_colidx, ᶜa_scalar_colidx, @@ -111,6 +116,7 @@ function edmfx_sgs_mass_flux_tendency!( (; ᶜρaʲs, ᶜρʲs, ᶠu³ʲs, ᶜKʲs, ᶜmseʲs, ᶜq_totʲs) = p.precomputed (; dt) = p ᶜJ = Fields.local_geometry_field(Y.c).J + ᶠJ = Fields.local_geometry_field(Y.f).J FT = eltype(Y) if p.atmos.edmfx_sgs_mass_flux @@ -140,6 +146,7 @@ function edmfx_sgs_mass_flux_tendency!( vertical_transport!( Yₜ.c.ρe_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³_diff_colidx, ᶜa_scalar_colidx, @@ -175,6 +182,7 @@ function edmfx_sgs_mass_flux_tendency!( vertical_transport!( Yₜ.c.ρq_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], ᶜρʲs.:($j)[colidx], ᶠu³_diff_colidx, ᶜa_scalar_colidx, diff --git a/src/prognostic_equations/implicit/implicit_solver.jl b/src/prognostic_equations/implicit/implicit_solver.jl index 40edcde99b..c32d0bcbd3 100644 --- a/src/prognostic_equations/implicit/implicit_solver.jl +++ b/src/prognostic_equations/implicit/implicit_solver.jl @@ -499,6 +499,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ᶜuₕ = Y.c.uₕ ᶠu₃ = Y.f.u₃ ᶜJ = Fields.local_geometry_field(Y.c).J + ᶠJ = Fields.local_geometry_field(Y.f).J ᶜgⁱʲ = Fields.local_geometry_field(Y.c).gⁱʲ ᶠgⁱʲ = Fields.local_geometry_field(Y.f).gⁱʲ @@ -507,14 +508,8 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) TD.gas_constant_air(thermo_params, ᶜts[colidx]) / TD.cv_m(thermo_params, ᶜts[colidx]) - if use_derivative(topography_flag) - @. ∂ᶜK_∂ᶜuₕ[colidx] = DiagonalMatrixRow( - adjoint(CTh(ᶜuₕ[colidx])) + - adjoint(ᶜinterp(ᶠu₃[colidx])) * g³ʰ(ᶜgⁱʲ[colidx]), - ) - else - @. ∂ᶜK_∂ᶜuₕ[colidx] = DiagonalMatrixRow(adjoint(CTh(ᶜuₕ[colidx]))) - end + @. ∂ᶜK_∂ᶜuₕ[colidx] = + DiagonalMatrixRow(adjoint(CTh(ᶜuₕ[colidx]) + CTh(ᶜinterp(ᶠu₃[colidx])))) @. ∂ᶜK_∂ᶠu₃[colidx] = ᶜinterp_matrix() ⋅ DiagonalMatrixRow(adjoint(CT3(ᶠu₃[colidx]))) + DiagonalMatrixRow(adjoint(CT3(ᶜuₕ[colidx]))) ⋅ ᶜinterp_matrix() @@ -524,7 +519,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) @. ᶜadvection_matrix[colidx] = -(ᶜadvdivᵥ_matrix()) ⋅ - DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρ[colidx])) + DiagonalMatrixRow(ᶠinterp(ᶜJ[colidx] * ᶜρ[colidx]) / ᶠJ[colidx]) if use_derivative(topography_flag) ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)] @@ -716,7 +711,7 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ᶠtmp = p.ᶠtemp_CT3 @. ᶠtmp[colidx] = CT3(unit_basis_vector_data(CT3, ᶠlg[colidx])) * - ᶠwinterp(ᶜJ[colidx], ᶜρ[colidx]) + ᶠinterp(ᶜJ[colidx] * ᶜρ[colidx]) / ᶠJ[colidx] @. ∂ᶜρqₚ_err_∂ᶜρqₚ[colidx] += dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠtmp) ⋅ ᶠright_bias_matrix() ⋅ @@ -802,14 +797,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ᶜρʲs.:(1)[colidx], ), ), - ), - ) ⋅ ᶠwinterp_matrix(ᶜJ[colidx]) ⋅ DiagonalMatrixRow( - ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / + ) / ᶠJ[colidx], + ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( + ᶜJ[colidx] * ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]) * e_int_v0, ) @. ᶠbidiagonal_matrix_ct3_2[colidx] = - DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅ + DiagonalMatrixRow( + ᶠinterp(ᶜJ[colidx] * ᶜρʲs.:(1)[colidx]) / ᶠJ[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, @@ -832,14 +828,15 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ᶜρʲs.:(1)[colidx], ), ), - ), - ) ⋅ ᶠwinterp_matrix(ᶜJ[colidx]) ⋅ DiagonalMatrixRow( - ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / + ) / ᶠJ[colidx], + ) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow( + ᶜJ[colidx] * ᶜkappa_mʲ[colidx] * (ᶜρʲs.:(1)[colidx])^2 / ((ᶜkappa_mʲ[colidx] + 1) * ᶜp[colidx]), ) @. ᶠbidiagonal_matrix_ct3_2[colidx] = - DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) ⋅ - ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅ + DiagonalMatrixRow( + ᶠinterp(ᶜJ[colidx] * ᶜρʲs.:(1)[colidx]) / ᶠJ[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]), @@ -852,8 +849,9 @@ function update_implicit_equation_jacobian!(A, Y, p, dtγ, colidx) ∂ᶜρaʲ_err_∂ᶜρaʲ = matrix[@name(c.sgsʲs.:(1).ρa), @name(c.sgsʲs.:(1).ρa)] @. ᶜadvection_matrix[colidx] = - -(ᶜadvdivᵥ_matrix()) ⋅ - DiagonalMatrixRow(ᶠwinterp(ᶜJ[colidx], ᶜρʲs.:(1)[colidx])) + -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow( + ᶠinterp(ᶜJ[colidx] * ᶜρʲs.:(1)[colidx]) / ᶠJ[colidx], + ) @. ∂ᶜρaʲ_err_∂ᶜρaʲ[colidx] = dtγ * ᶜadvection_matrix[colidx] ⋅ ᶠset_upwind_matrix_bcs(ᶠupwind_matrix(ᶠu³ʲs.:(1)[colidx])) ⋅ diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index e3b217dc5a..aeed0f2cfd 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -56,49 +56,64 @@ end # the implicit tendency function. Since dt >= dtγ, we can safely use dt for now. # TODO: Can we rewrite ᶠfct_boris_book and ᶠfct_zalesak so that their broadcast # expressions are less convoluted? -vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val, ᶜdivᵥ) = - vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜdivᵥ) -vertical_transport!(ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val) = - vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ) +vertical_transport!(ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val, ᶜdivᵥ) = + vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜdivᵥ) +vertical_transport!(ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding::Val) = + vertical_transport!(1, ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ) vertical_transport!( coeff::Int, ᶜρχₜ, ᶜJ, + ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt::Real, upwinding::Val, -) = vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, upwinding, ᶜadvdivᵥ) +) = vertical_transport!( + coeff, + ᶜρχₜ, + ᶜJ, + ᶠJ, + ᶜρ, + ᶠu³, + ᶜχ, + dt, + upwinding, + ᶜadvdivᵥ, +) -vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}, ᶜdivᵥ) = - @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ))) +vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}, ᶜdivᵥ) = + @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠu³ * ᶠinterp(ᶜχ) / ᶠJ)) vertical_transport!( coeff, ᶜρχₜ, ᶜJ, + ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order}, ᶜdivᵥ, -) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ))) +) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ) / ᶠJ)) vertical_transport!( coeff, ᶜρχₜ, ᶜJ, + ᶠJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order}, ᶜdivᵥ, -) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ))) +) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind3(ᶠu³, ᶜχ) / ᶠJ)) vertical_transport!( coeff, ᶜρχₜ, ᶜJ, + ᶠJ, ᶜρ, ᶠu³, ᶜχ, @@ -107,24 +122,34 @@ vertical_transport!( ᶜdivᵥ, ) = @. ᶜρχₜ += -coeff * (ᶜdivᵥ( - ᶠwinterp(ᶜJ, ᶜρ) * ( + ᶠinterp(ᶜJ * ᶜρ) * ( ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_boris_book( ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ), - ᶜχ / dt - ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ, + ᶜχ / dt - ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ) / ᶠJ) / ᶜρ, + ) + ) / ᶠJ, + )) +vertical_transport!( + coeff, + ᶜρχₜ, + ᶜJ, + ᶠJ, + ᶜρ, + ᶠu³, + ᶜχ, + dt, + ::Val{:zalesak}, + ᶜdivᵥ, +) = @. ᶜρχₜ += + -coeff * (ᶜdivᵥ( + ᶠinterp(ᶜJ * ᶜρ) * ( + ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_zalesak( + ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ), + ᶜχ / dt, + ᶜχ / dt - ᶜdivᵥ(ᶠinterp(ᶜJ * ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ) / ᶠJ) / ᶜρ, ) - ), + ) / ᶠJ, )) -vertical_transport!(coeff, ᶜρχₜ, ᶜJ, ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:zalesak}, ᶜdivᵥ) = - @. ᶜρχₜ += - -coeff * (ᶜdivᵥ( - ᶠwinterp(ᶜJ, ᶜρ) * ( - ᶠupwind1(ᶠu³, ᶜχ) + ᶠfct_zalesak( - ᶠupwind3(ᶠu³, ᶜχ) - ᶠupwind1(ᶠu³, ᶜχ), - ᶜχ / dt, - ᶜχ / dt - ᶜdivᵥ(ᶠwinterp(ᶜJ, ᶜρ) * ᶠupwind1(ᶠu³, ᶜχ)) / ᶜρ, - ) - ), - )) vertical_advection!(ᶜρχₜ, ᶠu³, ᶜχ, ::Val{:none}) = @. ᶜρχₜ -= ᶜadvdivᵥ(ᶠu³ * ᶠinterp(ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³) @@ -138,16 +163,18 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) (; dt) = p n = n_mass_flux_subdomains(turbconv_model) ᶜJ = Fields.local_geometry_field(Y.c).J + ᶠJ = Fields.local_geometry_field(Y.f).J (; ᶠgradᵥ_ᶜΦ, ᶜρ_ref, ᶜp_ref) = p.core (; ᶜh_tot, ᶜspecific, ᶠu³, ᶜp) = p.precomputed @. Yₜ.c.ρ[colidx] -= - ᶜdivᵥ(ᶠwinterp(ᶜJ[colidx], Y.c.ρ[colidx]) * ᶠu³[colidx]) + ᶜdivᵥ(ᶠinterp(ᶜJ[colidx] * Y.c.ρ[colidx]) * ᶠu³[colidx] / ᶠJ[colidx]) # Central advection of active tracers (e_tot and q_tot) vertical_transport!( Yₜ.c.ρe_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], Y.c.ρ[colidx], ᶠu³[colidx], ᶜh_tot[colidx], @@ -158,6 +185,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) vertical_transport!( Yₜ.c.ρq_tot[colidx], ᶜJ[colidx], + ᶠJ[colidx], Y.c.ρ[colidx], ᶠu³[colidx], ᶜspecific.q_tot[colidx], @@ -175,13 +203,13 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t, colidx) ᶠlg = Fields.local_geometry_field(Y.f) @. Yₜ.c.ρq_rai[colidx] -= ᶜprecipdivᵥ( CT3(unit_basis_vector_data(CT3, ᶠlg[colidx])) * - ᶠwinterp(ᶜJ[colidx], Y.c.ρ[colidx]) * - ᶠright_bias(-p.precomputed.ᶜwᵣ[colidx] * ᶜspecific.q_rai[colidx]), + ᶠinterp(ᶜJ[colidx] * Y.c.ρ[colidx]) * + ᶠright_bias(-p.precomputed.ᶜwᵣ[colidx] * ᶜspecific.q_rai[colidx]) / ᶠJ[colidx], ) @. Yₜ.c.ρq_sno[colidx] -= ᶜprecipdivᵥ( CT3(unit_basis_vector_data(CT3, ᶠlg[colidx])) * - ᶠwinterp(ᶜJ[colidx], Y.c.ρ[colidx]) * - ᶠright_bias(-p.precomputed.ᶜwₛ[colidx] * ᶜspecific.q_sno[colidx]), + ᶠinterp(ᶜJ[colidx] * Y.c.ρ[colidx]) * + ᶠright_bias(-p.precomputed.ᶜwₛ[colidx] * ᶜspecific.q_sno[colidx]) / ᶠJ[colidx], ) end diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl index 2a37106d55..5154caf855 100644 --- a/src/utils/utilities.jl +++ b/src/utils/utilities.jl @@ -61,11 +61,11 @@ function compute_kinetic!(κ::Fields.Field, uₕ::Fields.Field, uᵥ::Fields.Fie @assert eltype(uₕ) <: Union{C1, C2, C12} @assert eltype(uᵥ) <: C3 @. κ = - 1 / 2 * ( - dot(C123(uₕ), CT123(uₕ)) + - ᶜinterp(dot(C123(uᵥ), CT123(uᵥ))) + - 2 * dot(CT123(uₕ), ᶜinterp(C123(uᵥ))) - ) + ( + dot(uₕ, CT12(uₕ) + CT12(ᶜinterp(uᵥ))) + + dot(ᶜinterp(uᵥ), CT3(uₕ)) + + ᶜinterp(dot(uᵥ, CT3(uᵥ))) + ) / 2 end """