solvers_parabolic.jl

"""
  ParabolicFormulationBassiRebay1()

The classical BR1 flux from

- F. Bassi, S. Rebay (1997)
  A High-Order Accurate Discontinuous Finite Element Method for
  the Numerical Solution of the Compressible Navier-Stokes Equations
  [DOI: 10.1006/jcph.1996.5572](https://doi.org/10.1006/jcph.1996.5572)

A more detailed study of the BR1 scheme for the DGSEM can be found in
- G. J. Gassner, A. R. Winters, F. J. Hindenlang, D. Kopriva (2018)
  The BR1 Scheme is Stable for the Compressible Navier-Stokes Equations
  [DOI: 10.1007/s10915-018-0702-1](https://doi.org/10.1007/s10915-018-0702-1)

The BR1 scheme works well for convection-dominated problems, but may cause instabilities or 
reduced convergence for diffusion-dominated problems. 
In the latter case, the [`ParabolicFormulationLocalDG`](@ref) scheme is recommended.
"""
struct ParabolicFormulationBassiRebay1 end

"""
    flux_parabolic(u_ll, u_rr,
                   gradient_or_divergence, equations_parabolic,
                   parabolic_scheme::ParabolicFormulationBassiRebay1)

    flux_parabolic(u_ll, u_rr, normal_direction::AbstractVector,
                   gradient_or_divergence, equations_parabolic,
                   parabolic_scheme::ParabolicFormulationBassiRebay1)

This computes the classical BR1 flux. Since the interface flux for both the 
DG gradient and DG divergence under BR1 are identical, this function does 
not need to be specialized for `Gradient` and `Divergence`.

`normal_direction` is not used in the BR1 flux,
but is included as an argument for consistency with the [`ParabolicFormulationLocalDG`](@ref) flux,
which does use the `normal_direction` to compute the LDG "switch" on the generally non-Cartesian [`P4estMesh`](@ref).
"""
function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh`
                        gradient_or_divergence, equations_parabolic,
                        parabolic_scheme::ParabolicFormulationBassiRebay1)
    return 0.5f0 * (u_ll + u_rr)
end
# Version for `P4estMesh`
function flux_parabolic(u_ll, u_rr, normal_direction::AbstractVector,
                        gradient_or_divergence, equations_parabolic,
                        parabolic_scheme::ParabolicFormulationBassiRebay1)
    return 0.5f0 * (u_ll + u_rr)
end

"""
    ParabolicFormulationLocalDG(penalty_parameter)

The local DG (LDG) flux from "The Local Discontinuous Galerkin Method for Time-Dependent
Convection-Diffusion Systems" by Cockburn and Shu (1998).

The parabolic "upwinding" vector is currently implemented for `TreeMesh`; for all other mesh types,
the LDG solver is equivalent to [`ParabolicFormulationBassiRebay1`](@ref) with an LDG-type penalization.

- Cockburn and Shu (1998).
  The Local Discontinuous Galerkin Method for Time-Dependent
  Convection-Diffusion Systems
  [DOI: 10.1137/S0036142997316712](https://doi.org/10.1137/S0036142997316712)
"""
struct ParabolicFormulationLocalDG{P}
    penalty_parameter::P
end

"""
    ParabolicFormulationLocalDG()

The minimum dissipation local DG (LDG) flux from "An Analysis of the Minimal Dissipation Local 
Discontinuous Galerkin Method for Convection–Diffusion Problems" by Cockburn and Dong (2007). 
This scheme corresponds to an LDG parabolic "upwinding/downwinding" but no LDG penalty parameter. 
Cockburn and Dong proved that this scheme is still stable despite the zero penalty parameter. 

- Cockburn and Dong (2007)  
  An Analysis of the Minimal Dissipation Local Discontinuous 
  Galerkin Method for Convection–Diffusion Problems.
  [DOI: 10.1007/s10915-007-9130-3](https://doi.org/10.1007/s10915-007-9130-3)
"""
ParabolicFormulationLocalDG() = ParabolicFormulationLocalDG(nothing)

@doc raw"""
    flux_parabolic(u_ll, u_rr,
                   ::Gradient, equations_parabolic,
                   parabolic_scheme::ParabolicFormulationLocalDG)

    flux_parabolic(u_ll, u_rr, normal_direction,
                   ::Gradient, equations_parabolic,
                   parabolic_scheme::ParabolicFormulationLocalDG)

These fluxes computes the gradient and divergence interface fluxes for the 
local DG method. The local DG method uses an "upwind/downwind" flux for the 
gradient and divergence (i.e., if the gradient is upwinded, the divergence
must be downwinded in order to preserve symmetry and positive definiteness).
Here, we use the convention that the gradient flux is upwinded, thus we have
```math
f_{\text{gradient}} = u_{L}
```
on the Cartesian [`TreeMesh`](@ref).

For the [`P4estMesh`](@ref), the `normal_direction` is used to compute the LDG "switch" ``\sigma`` for the upwinding.
This is realized by selecting the sign of the maximum (in absolute value sense) normal direction component,
which corresponds to the "dominant" direction of the interface normal.
```math
i = \text{argmax} \{ \begin{pmatrix} \vert n_1 \vert \\ \vert n_2 \vert \\ \dots \end{pmatrix} \}
\sigma = \text{sign} (n_i)
```
"""
function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh`
                        ::Gradient, equations_parabolic,
                        parabolic_scheme::ParabolicFormulationLocalDG)
    # The LDG flux is {{f}} + beta * [[f]], where beta is the LDG "switch", 
    # which we set to -1 on the left and +1 on the right in 1D. The sign of the 
    # jump term should be opposite that of the sign used in the divergence flux. 
    # This is equivalent to setting the flux equal to `u_ll` for the gradient,
    # and `u_rr` for the divergence. 
    return u_ll # Use the upwind value for the gradient interface flux
end
# Version for `P4estMesh`
function flux_parabolic(u_ll, u_rr, normal_direction,
                        ::Gradient, equations_parabolic,
                        parabolic_scheme::ParabolicFormulationLocalDG)
    # Use "Upwind in dominant direction" for LDG switch
    abs_max_dir = argmax(abs.(normal_direction))
    ldg_switch = sign(normal_direction[abs_max_dir])
    return 0.5f0 * (u_ll + u_rr - ldg_switch * (u_rr - u_ll))
end

@doc raw"""
    flux_parabolic(u_ll, u_rr,
                   ::Divergence, equations_parabolic,
                   parabolic_scheme::ParabolicFormulationLocalDG)

    flux_parabolic(u_ll, u_rr, normal_direction,
                   ::Divergence, equations_parabolic,  
                   parabolic_scheme::ParabolicFormulationLocalDG)

These fluxes computes the gradient and divergence interface fluxes for the 
local DG method. The local DG method uses an "upwind/downwind" flux for the 
gradient and divergence (i.e., if the gradient is upwinded, the divergence
must be downwinded in order to preserve symmetry and positive definiteness).
Here, we use the convention that, because the gradient flux is upwinded, the divergence flux is downwinded.
Thus we have
```math
f_{\text{divergence}} = u_{R}
```
on the Cartesian [`TreeMesh`](@ref).

For the [`P4estMesh`](@ref), the `normal_direction` is used to compute the LDG "switch" ``\sigma`` for the downwinding.
This is realized by selecting the sign of the maximum (in absolute value sense) normal direction component,
which corresponds to the "dominant" direction of the interface normal.
```math
i = \text{argmax} \{ \begin{pmatrix} \vert n_1 \vert \\ \vert n_2 \vert \\ \dots \end{pmatrix} \}
\sigma = -\text{sign} (n_i)
```
"""
function flux_parabolic(u_ll, u_rr, # Version for `TreeMesh`
                        ::Divergence, equations_parabolic,
                        parabolic_scheme::ParabolicFormulationLocalDG)
    return u_rr # Use the downwind value for the divergence interface flux
end
# Version for `P4estMesh`
function flux_parabolic(u_ll, u_rr, normal_direction,
                        ::Divergence, equations_parabolic,
                        parabolic_scheme::ParabolicFormulationLocalDG)
    # Use "Downwind in dominant direction" for LDG switch
    abs_max_dir = argmax(abs.(normal_direction))
    ldg_switch = -sign(normal_direction[abs_max_dir])
    return 0.5f0 * (u_ll + u_rr - ldg_switch * (u_rr - u_ll))
end

default_parabolic_solver() = ParabolicFormulationBassiRebay1()