compressible_navier_stokes_3d.jl

# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin
#! format: noindent

@doc raw"""
    CompressibleNavierStokesDiffusion3D(equations; mu, Pr,
                                        gradient_variables=GradientVariablesPrimitive())

Contains the diffusion (i.e. parabolic) terms applied
to mass, momenta, and total energy together with the advective terms from
the [`CompressibleEulerEquations3D`](@ref).

- `equations`: instance of the [`CompressibleEulerEquations3D`](@ref)
- `mu`: dynamic viscosity,
- `Pr`: Prandtl number,
- `gradient_variables`: which variables the gradients are taken with respect to.
                        Defaults to [`GradientVariablesPrimitive()`](@ref).
                        For an entropy stable formulation, use [`GradientVariablesEntropy()`](@ref).

Fluid properties such as the dynamic viscosity ``\mu`` can be provided in any consistent unit system, e.g.,
[``\mu``] = kg m⁻¹ s⁻¹.
The viscosity ``\mu`` may be a constant or a function of the current state, e.g.,
depending on temperature (Sutherland's law): ``\mu = \mu(T)``.
In the latter case, the function `mu` needs to have the signature `mu(u, equations)`.

The particular form of the compressible Navier-Stokes implemented is
```math
\frac{\partial}{\partial t}
\begin{pmatrix}
\rho \\ \rho \mathbf{v} \\ \rho e_{\text{total}}
\end{pmatrix}
+
\nabla \cdot
\begin{pmatrix}
 \rho \mathbf{v} \\ \rho \mathbf{v}\mathbf{v}^T + p \underline{I} \\ (\rho e_{\text{total}} + p) \mathbf{v}
\end{pmatrix}
=
\nabla \cdot
\begin{pmatrix}
0 \\ \underline{\tau} \\ \underline{\tau}\mathbf{v} - \mathbf{q}
\end{pmatrix}
```
where the system is closed with the ideal gas assumption giving
```math
p = (\gamma - 1) \left( \rho e_{\text{total}} - \frac{1}{2} \rho (v_1^2+v_2^2+v_3^2) \right)
```
as the pressure. The value of the adiabatic constant `gamma` is taken from the [`CompressibleEulerEquations3D`](@ref).
The terms on the right hand side of the system above
are built from the viscous stress tensor
```math
\underline{\tau} = \mu \left(\nabla\mathbf{v} + \left(\nabla\mathbf{v}\right)^T\right) - \frac{2}{3} \mu \left(\nabla\cdot\mathbf{v}\right)\underline{I}
```
where ``\underline{I}`` is the ``3\times 3`` identity matrix and the heat flux is
```math
\mathbf{q} = -\kappa\nabla\left(T\right),\quad T = \frac{p}{R\rho}
```
where ``T`` is the temperature and ``\kappa`` is the thermal conductivity for Fourier's law.
Under the assumption that the gas has a constant Prandtl number,
the thermal conductivity is
```math
\kappa = \frac{\gamma \mu R}{(\gamma - 1)\textrm{Pr}}.
```
From this combination of temperature ``T`` and thermal conductivity ``\kappa`` we see
that the gas constant `R` cancels and the heat flux becomes
```math
\mathbf{q} = -\kappa\nabla\left(T\right) = -\frac{\gamma \mu}{(\gamma - 1)\textrm{Pr}}\nabla\left(\frac{p}{\rho}\right)
```
which is the form implemented below in the [`flux`](@ref) function.

In three spatial dimensions we require gradients for four quantities, e.g.,
primitive quantities
```math
\nabla v_1,\, \nabla v_2,\, \nabla v_3,\, \nabla T
```
or the entropy variables
```math
\nabla w_2,\, \nabla w_3,\, \nabla w_4\, \nabla w_5
```
where
```math
w_2 = \frac{\rho v_1}{p},\, w_3 = \frac{\rho v_2}{p},\, w_4 = \frac{\rho v_3}{p},\, w_5 = -\frac{\rho}{p}
```
"""
struct CompressibleNavierStokesDiffusion3D{GradientVariables, RealT <: Real, Mu,
                                           E <: AbstractCompressibleEulerEquations{3}} <:
       AbstractCompressibleNavierStokesDiffusion{3, 5, GradientVariables}
    # TODO: parabolic
    # Add NGRADS as a type parameter here and in AbstractEquationsParabolic, add `ngradients(...)` accessor function

    mu::Mu                     # viscosity
    Pr::RealT                  # Prandtl number
    kappa::RealT               # thermal diffusivity for Fourier's law
    max_4over3_kappa::RealT    # max(4/3, kappa) used for parabolic cfl => `max_diffusivity`

    equations_hyperbolic::E    # CompressibleEulerEquations3D
    gradient_variables::GradientVariables # GradientVariablesPrimitive or GradientVariablesEntropy
end

# default to primitive gradient variables
function CompressibleNavierStokesDiffusion3D(equations::CompressibleEulerEquations3D;
                                             mu, Prandtl,
                                             gradient_variables = GradientVariablesPrimitive())
    @unpack gamma, inv_gamma_minus_one = equations

    Pr = promote_type(typeof(gamma), typeof(Prandtl))(Prandtl)
    # Under the assumption of constant Prandtl number the thermal conductivity
    # constant is kappa = gamma μ / ((gamma-1) Prandtl).
    # Important note! Factor of μ is accounted for later in `flux`.
    # This avoids recomputation of kappa for non-constant μ.
    kappa = gamma * inv_gamma_minus_one / Pr

    return CompressibleNavierStokesDiffusion3D{typeof(gradient_variables),
                                               typeof(Pr), typeof(mu),
                                               typeof(equations)}(mu, Pr, kappa,
                                                                  max(4 / 3, kappa),
                                                                  equations,
                                                                  gradient_variables)
end

# TODO: parabolic
# This is the flexibility a user should have to select the different gradient variable types
# varnames(::typeof(cons2prim)   , ::CompressibleNavierStokesDiffusion3D) = ("v1", "v2", "v3", "T")
# varnames(::typeof(cons2entropy), ::CompressibleNavierStokesDiffusion3D) = ("w2", "w3", "w4", "w5")

function varnames(variable_mapping,
                  equations_parabolic::CompressibleNavierStokesDiffusion3D)
    return varnames(variable_mapping, equations_parabolic.equations_hyperbolic)
end

# we specialize this function to compute gradients of primitive variables instead of
# conservative variables.
function gradient_variable_transformation(::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    return cons2prim
end
function gradient_variable_transformation(::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    return cons2entropy
end

# Explicit formulas for the diffusive Navier-Stokes fluxes are available, e.g., in Section 2
# of the paper by Rueda-Ramírez, Hennemann, Hindenlang, Winters, and Gassner
# "An Entropy Stable Nodal Discontinuous Galerkin Method for the resistive
#  MHD Equations. Part II: Subcell Finite Volume Shock Capturing"
# where one sets the magnetic field components equal to 0.
function flux(u, gradients, orientation::Integer,
              equations::CompressibleNavierStokesDiffusion3D)
    # Here, `u` is assumed to be the "transformed" variables specified by `gradient_variable_transformation`.
    _, v1, v2, v3, _ = convert_transformed_to_primitive(u, equations)
    # Here `gradients` is assumed to contain the gradients of the primitive variables (rho, v1, v2, v3, T)
    # either computed directly or reverse engineered from the gradient of the entropy variables
    # by way of the `convert_gradient_variables` function.
    _, dv1dx, dv2dx, dv3dx, dTdx = convert_derivative_to_primitive(u, gradients[1],
                                                                   equations)
    _, dv1dy, dv2dy, dv3dy, dTdy = convert_derivative_to_primitive(u, gradients[2],
                                                                   equations)
    _, dv1dz, dv2dz, dv3dz, dTdz = convert_derivative_to_primitive(u, gradients[3],
                                                                   equations)

    # Components of viscous stress tensor

    # Diagonal parts
    # (4 * (v1)_x / 3 - 2 * ((v2)_y + (v3)_z)) / 3)
    tau_11 = (4 * dv1dx - 2 * (dv2dy + dv3dz)) / 3
    # (4 * (v2)_y / 3 - 2 * ((v1)_x + (v3)_z) / 3)
    tau_22 = (4 * dv2dy - 2 * (dv1dx + dv3dz)) / 3
    # (4 * (v3)_z / 3 - 2 * ((v1)_x + (v2)_y) / 3)
    tau_33 = (4 * dv3dz - 2 * (dv1dx + dv2dy)) / 3

    # Off diagonal parts, exploit that stress tensor is symmetric
    # ((v1)_y + (v2)_x)
    tau_12 = dv1dy + dv2dx # = tau_21
    # ((v1)_z + (v3)_x)
    tau_13 = dv1dz + dv3dx # = tau_31
    # ((v2)_z + (v3)_y)
    tau_23 = dv2dz + dv3dy # = tau_32

    # Fourier's law q = -kappa * grad(T) = -kappa * grad(p / (R rho))
    # with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr)
    # Note, the gas constant cancels under this formulation, so it is not present
    # in the implementation
    q1 = equations.kappa * dTdx
    q2 = equations.kappa * dTdy
    q3 = equations.kappa * dTdz

    # In the simplest cases, the user passed in `mu` or `mu()`
    # (which returns just a constant) but
    # more complex functions like Sutherland's law are possible.
    # `dynamic_viscosity` is a helper function that handles both cases
    # by dispatching on the type of `equations.mu`.
    mu = dynamic_viscosity(u, equations)

    if orientation == 1
        # parabolic flux components in the x-direction
        f1 = 0
        f2 = tau_11 * mu
        f3 = tau_12 * mu
        f4 = tau_13 * mu
        f5 = (v1 * tau_11 + v2 * tau_12 + v3 * tau_13 + q1) * mu

        return SVector(f1, f2, f3, f4, f5)
    elseif orientation == 2
        # parabolic flux components in the y-direction
        # Note, symmetry is exploited for tau_12 = tau_21
        g1 = 0
        g2 = tau_12 * mu # tau_21 * mu
        g3 = tau_22 * mu
        g4 = tau_23 * mu
        g5 = (v1 * tau_12 + v2 * tau_22 + v3 * tau_23 + q2) * mu

        return SVector(g1, g2, g3, g4, g5)
    else # if orientation == 3
        # parabolic flux components in the z-direction
        # Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32
        h1 = 0
        h2 = tau_13 * mu # tau_31 * mu
        h3 = tau_23 * mu # tau_32 * mu
        h4 = tau_33 * mu
        h5 = (v1 * tau_13 + v2 * tau_23 + v3 * tau_33 + q3) * mu

        return SVector(h1, h2, h3, h4, h5)
    end
end

@doc raw"""
    max_diffusivity(u, equations_parabolic::CompressibleNavierStokesDiffusion3D)

# Returns
- `dynamic_viscosity(u, equations_parabolic) / u[1] * equations_parabolic.max_4over3_kappa`
where `max_4over3_kappa = max(4/3, kappa)` is computed in the constructor.

For the diffusive estimate we use the eigenvalues of the diffusivity matrix,
as suggested in Section 3.5 of
- Krais et. al (2021)
  FLEXI: A high order discontinuous Galerkin framework for hyperbolic–parabolic conservation laws
  [DOI: 10.1016/j.camwa.2020.05.004](https://doi.org/10.1016/j.camwa.2020.05.004)

For details on the derivation of eigenvalues of the diffusivity matrix
for the compressible Navier-Stokes equations see for instance
- Richard P. Dwight (2006)
  Efficiency improvements of RANS-based analysis and optimization using implicit and adjoint methods on unstructured grids
  PhD Thesis, University of Manchester
  https://elib.dlr.de/50794/1/rdwight-PhDThesis-ImplicitAndAdjoint.pdf
  See especially equations (2.79), (3.24), and (3.25) from Chapter 3.2.3

The eigenvalues of the diffusivity matrix in 3D are
``-\frac{\mu}{\rho} \{0, 4/3, 1, 1, \kappa\}``
and thus the largest absolute eigenvalue is
``\frac{\mu}{\rho} \max(4/3, \kappa)``.
"""
@inline function max_diffusivity(u,
                                 equations_parabolic::CompressibleNavierStokesDiffusion3D)
    # See for instance also the computation in FLUXO:
    # https://github.com/project-fluxo/fluxo/blob/c7e0cc9b7fd4569dcab67bbb6e5a25c0a84859f1/src/equation/navierstokes/calctimestep.f90#L122-L128
    #
    # Accordingly, the spectral radius/largest absolute eigenvalue can be computed as:
    return dynamic_viscosity(u, equations_parabolic) / u[1] *
           equations_parabolic.max_4over3_kappa
end

# Convert conservative variables to primitive
@inline function cons2prim(u, equations::CompressibleNavierStokesDiffusion3D)
    rho, rho_v1, rho_v2, rho_v3, _ = u

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    v3 = rho_v3 / rho
    T = temperature(u, equations)

    return SVector(rho, v1, v2, v3, T)
end

# Convert conservative variables to entropy
# TODO: parabolic. We can improve efficiency by not computing w_1, which involves logarithms
# This can be done by specializing `cons2entropy` and `entropy2cons` to `CompressibleNavierStokesDiffusion2D`,
# but this may be confusing to new users.
function cons2entropy(u, equations::CompressibleNavierStokesDiffusion3D)
    return cons2entropy(u, equations.equations_hyperbolic)
end
function entropy2cons(w, equations::CompressibleNavierStokesDiffusion3D)
    return entropy2cons(w, equations.equations_hyperbolic)
end

# the `flux` function takes in transformed variables `u` which depend on the type of the gradient variables.
# For CNS, it is simplest to formulate the parabolic terms in primitive variables, so we transform the transformed
# variables into primitive variables.
@inline function convert_transformed_to_primitive(u_transformed,
                                                  equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    return u_transformed
end

# TODO: parabolic. Make this more efficient!
@inline function convert_transformed_to_primitive(u_transformed,
                                                  equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    # note: this uses CompressibleNavierStokesDiffusion3D versions of cons2prim and entropy2cons
    return cons2prim(entropy2cons(u_transformed, equations), equations)
end

# Takes the solution values `u` and gradient of the entropy variables (w_2, w_3, w_4, w_5) and
# reverse engineers the gradients to be terms of the primitive variables (v1, v2, v3, T).
# Helpful because then the diffusive fluxes have the same form as on paper.
# Note, the first component of `gradient_entropy_vars` contains gradient(rho) which is unused.
# TODO: parabolic; entropy stable parabolic terms
@inline function convert_derivative_to_primitive(u, gradient,
                                                 ::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    return gradient
end

# the first argument is always the "transformed" variables.
@inline function convert_derivative_to_primitive(w, gradient_entropy_vars,
                                                 equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})

    # TODO: parabolic. This is inefficient to pass in transformed variables but then transform them back.
    # We can fix this if we directly compute v1, v2, v3, T from the entropy variables
    u = entropy2cons(w, equations) # calls a "modified" entropy2cons defined for CompressibleNavierStokesDiffusion3D
    rho, rho_v1, rho_v2, rho_v3, _ = u

    v1 = rho_v1 / rho
    v2 = rho_v2 / rho
    v3 = rho_v3 / rho
    T = temperature(u, equations)

    return SVector(gradient_entropy_vars[1],
                   T * (gradient_entropy_vars[2] + v1 * gradient_entropy_vars[5]), # grad(u) = T*(grad(w_2)+v1*grad(w_5))
                   T * (gradient_entropy_vars[3] + v2 * gradient_entropy_vars[5]), # grad(v) = T*(grad(w_3)+v2*grad(w_5))
                   T * (gradient_entropy_vars[4] + v3 * gradient_entropy_vars[5]), # grad(v) = T*(grad(w_4)+v3*grad(w_5))
                   T * T * gradient_entropy_vars[5])
end

# This routine is required because `prim2cons` is called in `initial_condition`, which
# is called with `equations::CompressibleEulerEquations3D`. This means it is inconsistent
# with `cons2prim(..., ::CompressibleNavierStokesDiffusion3D)` as defined above.
# TODO: parabolic. Is there a way to clean this up?
@inline function prim2cons(u, equations::CompressibleNavierStokesDiffusion3D)
    return prim2cons(u, equations.equations_hyperbolic)
end

"""
    temperature(u, equations::CompressibleNavierStokesDiffusion3D)

Compute the temperature from the conservative variables `u`.
In particular, this assumes a specific gas constant ``R = 1``:
```math
T = \\frac{p}{\\rho}
```
"""
@inline function temperature(u, equations::CompressibleNavierStokesDiffusion3D)
    rho, rho_v1, rho_v2, rho_v3, rho_e_total = u
    @unpack gamma = equations

    p = (gamma - 1) *
        (rho_e_total - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho)
    T = p / rho # Corresponds to a specific gas constant R = 1
    return T
end

@inline function velocity(u, equations::CompressibleNavierStokesDiffusion3D)
    rho = u[1]
    v1 = u[2] / rho
    v2 = u[3] / rho
    v3 = u[4] / rho
    return SVector(v1, v2, v3)
end

@doc raw"""
    enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D)

Computes the (node-wise) enstrophy, defined as
```math
    \mathcal{E} = \frac{1}{2} \rho \boldsymbol{\omega} \cdot \boldsymbol{\omega}
```
where ``\boldsymbol{\omega} = \nabla \times \boldsymbol{v}`` is the [`vorticity`](@ref).
In 3D, ``\boldsymbol{\omega}`` is a full three-component vector.
"""
@inline function enstrophy(u, gradients, equations::CompressibleNavierStokesDiffusion3D)
    # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v

    omega = vorticity(u, gradients, equations)
    return 0.5f0 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2)
end

@doc raw"""
    vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D)

Computes the (node-wise) vorticity, defined in 3D as
```math
    \boldsymbol{\omega} = \nabla \times \boldsymbol{v} =
    \begin{pmatrix}
        \frac{\partial v_3}{\partial y} - \frac{\partial v_2}{\partial z} \\
        \frac{\partial v_1}{\partial z} - \frac{\partial v_3}{\partial x} \\
        \frac{\partial v_2}{\partial x} - \frac{\partial v_1}{\partial y}
    \end{pmatrix}
```
"""
@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D)
    # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function.
    _, _, dv2dx, dv3dx, _ = convert_derivative_to_primitive(u, gradients[1], equations)
    _, dv1dy, _, dv3dy, _ = convert_derivative_to_primitive(u, gradients[2], equations)
    _, dv1dz, dv2dz, _, _ = convert_derivative_to_primitive(u, gradients[3], equations)

    return SVector(dv3dy - dv2dz, dv1dz - dv3dx, dv2dx - dv1dy)
end

@inline function vorticity(u, gradients,
                           equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    # Need to convert to entropy variables first for `convert_derivative_to_primitive` to work correctly.
    w = cons2entropy(u, equations)

    # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function.
    _, _, dv2dx, dv3dx, _ = convert_derivative_to_primitive(w, gradients[1], equations)
    _, dv1dy, _, dv3dy, _ = convert_derivative_to_primitive(w, gradients[2], equations)
    _, dv1dz, dv2dz, _, _ = convert_derivative_to_primitive(w, gradients[3], equations)

    return SVector(dv3dy - dv2dz, dv1dz - dv3dx, dv2dx - dv1dy)
end

@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Adiabatic})(flux_inner,
                                                                                      u_inner,
                                                                                      normal::AbstractVector,
                                                                                      x,
                                                                                      t,
                                                                                      operator_type::Gradient,
                                                                                      equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x,
                                                                                        t,
                                                                                        equations)
    return SVector(u_inner[1], v1, v2, v3, u_inner[5])
end

@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Adiabatic})(flux_inner,
                                                                                      u_inner,
                                                                                      normal::AbstractVector,
                                                                                      x,
                                                                                      t,
                                                                                      operator_type::Divergence,
                                                                                      equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x,
                                                                                                           t,
                                                                                                           equations)
    v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x,
                                                                                        t,
                                                                                        equations)
    _, tau_1n, tau_2n, tau_3n, _ = flux_inner # extract fluxes for 2nd, 3rd, and 4th equations
    normal_energy_flux = v1 * tau_1n + v2 * tau_2n + v3 * tau_3n + normal_heat_flux
    return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4],
                   normal_energy_flux)
end

@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Isothermal})(flux_inner,
                                                                                       u_inner,
                                                                                       normal::AbstractVector,
                                                                                       x,
                                                                                       t,
                                                                                       operator_type::Gradient,
                                                                                       equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x,
                                                                                        t,
                                                                                        equations)
    T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t,
                                                                                equations)
    return SVector(u_inner[1], v1, v2, v3, T)
end

@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Isothermal})(flux_inner,
                                                                                       u_inner,
                                                                                       normal::AbstractVector,
                                                                                       x,
                                                                                       t,
                                                                                       operator_type::Divergence,
                                                                                       equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    return flux_inner
end

# specialized BC impositions for GradientVariablesEntropy.

# This should return a SVector containing the boundary values of entropy variables.
# Here, `w_inner` are the transformed variables (e.g., entropy variables).
#
# Taken from "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions
#             for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022.
# DOI: 10.1016/j.jcp.2021.110723
@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Adiabatic})(flux_inner,
                                                                                      w_inner,
                                                                                      normal::AbstractVector,
                                                                                      x,
                                                                                      t,
                                                                                      operator_type::Gradient,
                                                                                      equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x,
                                                                                        t,
                                                                                        equations)
    negative_rho_inv_p = w_inner[5] # w_5 = -rho / p
    return SVector(w_inner[1], -v1 * negative_rho_inv_p, -v2 * negative_rho_inv_p,
                   -v3 * negative_rho_inv_p, negative_rho_inv_p)
end

# this is actually identical to the specialization for GradientVariablesPrimitive, but included for completeness.
@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Adiabatic})(flux_inner,
                                                                                      w_inner,
                                                                                      normal::AbstractVector,
                                                                                      x,
                                                                                      t,
                                                                                      operator_type::Divergence,
                                                                                      equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x,
                                                                                                           t,
                                                                                                           equations)
    v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x,
                                                                                        t,
                                                                                        equations)
    _, tau_1n, tau_2n, tau_3n, _ = flux_inner # extract fluxes for 2nd, 3rd, and 4th equations
    normal_energy_flux = v1 * tau_1n + v2 * tau_2n + v3 * tau_3n + normal_heat_flux
    return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4],
                   normal_energy_flux)
end

@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Isothermal})(flux_inner,
                                                                                       w_inner,
                                                                                       normal::AbstractVector,
                                                                                       x,
                                                                                       t,
                                                                                       operator_type::Gradient,
                                                                                       equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    v1, v2, v3 = boundary_condition.boundary_condition_velocity.boundary_value_function(x,
                                                                                        t,
                                                                                        equations)
    T = boundary_condition.boundary_condition_heat_flux.boundary_value_function(x, t,
                                                                                equations)

    # the entropy variables w2 = rho * v1 / p = v1 / T = -v1 * w5. Similarly for w3 and w4
    w5 = -1 / T
    return SVector(w_inner[1], -v1 * w5, -v2 * w5, -v3 * w5, w5)
end

@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:NoSlip,
                                                                        <:Isothermal})(flux_inner,
                                                                                       w_inner,
                                                                                       normal::AbstractVector,
                                                                                       x,
                                                                                       t,
                                                                                       operator_type::Divergence,
                                                                                       equations::CompressibleNavierStokesDiffusion3D{GradientVariablesEntropy})
    return SVector(flux_inner[1], flux_inner[2], flux_inner[3], flux_inner[4],
                   flux_inner[5])
end

# Computes the mirror velocity across a symmetry plane which enforces
# a tangential velocity that is aligned with the symmetry plane, i.e.,
# which is normal to the `normal_direction`.
# See also `boundary_condition_slip_wall` and `rotate_to_x`.
@inline function velocity_symmetry_plane(normal_direction::AbstractVector, v1, v2, v3)
    norm_ = norm(normal_direction)
    normal = normal_direction / norm_

    # Compute alignment of velocity with normal direction
    v_normal = v1 * normal[1] + v2 * normal[2] + v3 * normal[3]

    v1_outer = v1 - 2 * v_normal * normal[1]
    v2_outer = v2 - 2 * v_normal * normal[2]
    v3_outer = v3 - 2 * v_normal * normal[3]

    return v1_outer, v2_outer, v3_outer
end

# Note: This should be used with `boundary_condition_slip_wall` for the hyperbolic (Euler) part.
@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:Slip,
                                                                        <:Adiabatic})(flux_inner,
                                                                                      u_inner,
                                                                                      normal::AbstractVector,
                                                                                      x,
                                                                                      t,
                                                                                      operator_type::Gradient,
                                                                                      equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    v1_outer, v2_outer, v3_outer = velocity_symmetry_plane(normal,
                                                           u_inner[2],
                                                           u_inner[3],
                                                           u_inner[4])

    return SVector(u_inner[1], v1_outer, v2_outer, v3_outer, u_inner[5])
end

# Note: This should be used with `boundary_condition_slip_wall` for the hyperbolic (Euler) part.
@inline function (boundary_condition::BoundaryConditionNavierStokesWall{<:Slip,
                                                                        <:Adiabatic})(flux_inner,
                                                                                      u_inner,
                                                                                      normal::AbstractVector,
                                                                                      x,
                                                                                      t,
                                                                                      operator_type::Divergence,
                                                                                      equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    normal_heat_flux = boundary_condition.boundary_condition_heat_flux.boundary_value_normal_flux_function(x,
                                                                                                           t,
                                                                                                           equations)
    # Normal stresses should be 0. This implies also that `normal_energy_flux = normal_heat_flux`.
    # For details, see Section 4.2 of
    # "Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions
    #  for the compressible Navier-Stokes equations" by Chan, Lin, Warburton 2022.
    # DOI: 10.1016/j.jcp.2021.110723
    return SVector(flux_inner[1], 0, 0, 0, normal_heat_flux)
end

# Dirichlet Boundary Condition for e.g. P4est mesh
@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
                                                                  u_inner,
                                                                  normal::AbstractVector,
                                                                  x, t,
                                                                  operator_type::Gradient,
                                                                  equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    # BCs are usually specified as conservative variables so we convert them to primitive variables
    #  because the gradients are assumed to be with respect to the primitive variables
    u_boundary = boundary_condition.boundary_value_function(x, t, equations)

    return cons2prim(u_boundary, equations)
end

@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner,
                                                                  u_inner,
                                                                  normal::AbstractVector,
                                                                  x, t,
                                                                  operator_type::Divergence,
                                                                  equations::CompressibleNavierStokesDiffusion3D{GradientVariablesPrimitive})
    # for Dirichlet boundary conditions, we do not impose any conditions on the parabolic fluxes
    return flux_inner
end
end # @muladd