Adding parabolic solver option

src/Trixi.jl

@@ -105,7 +105,7 @@ include("auxiliary/p4est.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
 include("solvers/solvers.jl")
-include("equations/parabolic_equations.jl") # these depend on parabolic solver types
+include("equations/equations_parabolic.jl") # these depend on parabolic solver types
 include("semidiscretization/semidiscretization.jl")
 include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
@@ -223,7 +223,7 @@ export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure
 export DGMulti, estimate_dt, DGMultiMesh, GaussSBP
 export VertexMappedMesh # TODO: DGMulti, v0.5. Remove deprecated VertexMappedMesh in next release
 
-export BR1, LDG
+export ViscousFormulationBassiRebay1, ViscousFormulationLocalDG
 
 # Visualization-related exports
 export PlotData1D, PlotData2D, ScalarPlotData2D, getmesh, adapt_to_mesh_level!, adapt_to_mesh_level

src/equations/laplace_diffusion_2d.jl

@@ -21,7 +21,10 @@ function flux(u, grad_u, equations::LaplaceDiffusion2D)
   return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy)
 end
 
-function penalty(u_outer, u_inner, inv_h, equations::LaplaceDiffusion2D, dg::LDG)
+# TODO: should this remain in the equations file, be moved to solvers, or live in the elixir?
+# The penalization depends on the solver, but also depends explicitly on physical parameters,
+# and would probably need to be specialized for every different equation.
+function penalty(u_outer, u_inner, inv_h, equations::LaplaceDiffusion2D, dg::ViscousFormulationLocalDG)
   return dg.penalty_parameter * (u_outer - u_inner) * equations.diffusivity * inv_h
 end
 

src/solvers/solvers.jl

@@ -7,7 +7,7 @@
 
 include("dg.jl")
 include("dgmulti.jl")
-include("parabolic_solvers.jl")
+include("solvers_parabolic.jl")
 
 
 end # @muladd

src/solvers/solvers_parabolic.jl

@@ -0,0 +1,38 @@
+"""
+  ViscousFormulationBassiRebay1()
+
+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)
+"""
+struct ViscousFormulationBassiRebay1 end
+
+# no penalization for a BR1 parabolic solver
+function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, boundary_conditions,
+                               mesh, equations::AbstractEquationsParabolic, dg::DGMulti,
+                               dg_parabolic::ViscousFormulationBassiRebay1, cache, cache_parabolic)
+  return nothing
+end
+
+"""
+    ViscousFormulationLocalDG(penalty_parameter)
+
+The local DG (LDG) flux from "The Local Discontinuous Galerkin Method for Time-Dependent
+Convection-Diffusion Systems" by Cockburn and Shu (1998).
+
+Note that, since this implementation does not involve the parabolic "upwinding" vector,
+the LDG solver is equivalent to [`ViscousFormulationBassiRebay1`](@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 ViscousFormulationLocalDG{P}
+  penalty_parameter::P
+end
+
+default_parabolic_solver() = ViscousFormulationBassiRebay1()