trixi-framework · jlchan · Jun 2, 2022 · May 31, 2022 · May 31, 2022 · May 31, 2022
diff --git a/src/Trixi.jl b/src/Trixi.jl
@@ -105,6 +105,7 @@ include("auxiliary/p4est.jl")
 include("equations/equations.jl")
 include("meshes/meshes.jl")
 include("solvers/solvers.jl")
+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")
@@ -222,6 +223,8 @@ 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 ViscousFormulationBassiRebay1, ViscousFormulationLocalDG
+
 # Visualization-related exports
 export PlotData1D, PlotData2D, ScalarPlotData2D, getmesh, adapt_to_mesh_level!, adapt_to_mesh_level
 

diff --git a/src/equations/equations.jl b/src/equations/equations.jl
@@ -370,8 +370,4 @@ include("acoustic_perturbation_2d.jl")
 
 abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
 
-# Linear scalar diffusion for use in linear scalar advection-diffusion problems
-abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end
-include("laplace_diffusion_2d.jl")
-
 end # @muladd
diff --git a/src/equations/equations_parabolic.jl b/src/equations/equations_parabolic.jl
@@ -0,0 +1,3 @@
+# Linear scalar diffusion for use in linear scalar advection-diffusion problems
+abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end
+include("laplace_diffusion_2d.jl")
diff --git a/src/equations/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl
@@ -21,6 +21,13 @@ function flux(u, grad_u, equations::LaplaceDiffusion2D)
   return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy)
 end
 
+# 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
+
 # Dirichlet-type boundary condition for use with a parabolic solver in weak form
 @inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector,
                                                                   x, t, operator_type::Gradient,

diff --git a/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl b/src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl
@@ -13,7 +13,7 @@ of a mixed hyperbolic-parabolic conservation law.
 """
 struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition,
                                              BoundaryConditions, BoundaryConditionsParabolic,
-                                             SourceTerms, Solver, Cache, CacheParabolic} <: AbstractSemidiscretization
+                                             SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic} <: AbstractSemidiscretization
 
   mesh::Mesh
 
@@ -30,17 +30,17 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic
   source_terms::SourceTerms
 
   solver::Solver
-  # TODO: introduce `solver_parabolic` for future specialization.
+  solver_parabolic::SolverParabolic
 
   cache::Cache
   cache_parabolic::CacheParabolic
 
   performance_counter::PerformanceCounter
 
-  function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}(
+  function SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}(
       mesh::Mesh, equations::Equations, equations_parabolic::EquationsParabolic, initial_condition::InitialCondition,
       boundary_conditions::BoundaryConditions, boundary_conditions_parabolic::BoundaryConditionsParabolic,
-      source_terms::SourceTerms, solver::Solver, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, Cache, CacheParabolic}
+      source_terms::SourceTerms, solver::Solver, solver_parabolic::SolverParabolic, cache::Cache, cache_parabolic::CacheParabolic) where {Mesh, Equations, EquationsParabolic, InitialCondition, BoundaryConditions, BoundaryConditionsParabolic, SourceTerms, Solver, SolverParabolic, Cache, CacheParabolic}
     @assert ndims(mesh) == ndims(equations)
 
     # Todo: assert nvariables(equations)==nvariables(equations_parabolic)
@@ -49,12 +49,13 @@ struct SemidiscretizationHyperbolicParabolic{Mesh, Equations, EquationsParabolic
 
     new(mesh, equations, equations_parabolic, initial_condition,
         boundary_conditions, boundary_conditions_parabolic,
-        source_terms, solver, cache, cache_parabolic, performance_counter)
+        source_terms, solver, solver_parabolic, cache, cache_parabolic, performance_counter)
   end
 end
 
 """
     SemidiscretizationHyperbolicParabolic(mesh, both_equations, initial_condition, solver;
+                                          solver_parabolic=default_parabolic_solver(),
                                           source_terms=nothing,
                                           both_boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic),
                                           RealT=real(solver),
@@ -65,6 +66,7 @@ Construct a semidiscretization of a hyperbolic-parabolic PDE.
 """
 function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple,
                                                initial_condition, solver;
+                                               solver_parabolic=default_parabolic_solver(),
                                                source_terms=nothing,
                                                boundary_conditions=(boundary_condition_periodic, boundary_condition_periodic),
                                                # `RealT` is used as real type for node locations etc.
@@ -77,7 +79,7 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations::Tuple,
   initial_hyperbolic_cache, initial_cache_parabolic = initial_caches
 
   return SemidiscretizationHyperbolicParabolic(mesh, equations_hyperbolic, equations_parabolic,
-                                               initial_condition, solver; source_terms,
+                                               initial_condition, solver; solver_parabolic, source_terms,
                                                boundary_conditions=boundary_conditions_hyperbolic,
                                                boundary_conditions_parabolic=boundary_conditions_parabolic,
                                                RealT, uEltype, initial_cache=initial_hyperbolic_cache,
@@ -86,6 +88,7 @@ end
 
 function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabolic,
                                                initial_condition, solver;
+                                               solver_parabolic=default_parabolic_solver(),
                                                source_terms=nothing,
                                                boundary_conditions=boundary_condition_periodic,
                                                boundary_conditions_parabolic=boundary_condition_periodic,
@@ -99,15 +102,15 @@ function SemidiscretizationHyperbolicParabolic(mesh, equations, equations_parabo
   _boundary_conditions = digest_boundary_conditions(boundary_conditions, mesh, solver, cache)
   _boundary_conditions_parabolic = digest_boundary_conditions(boundary_conditions_parabolic, mesh, solver, cache)
 
-  cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, RealT, uEltype)...,
+  cache_parabolic = (; create_cache_parabolic(mesh, equations_parabolic, solver, solver_parabolic, RealT, uEltype)...,
                                               initial_cache_parabolic...)
 
   SemidiscretizationHyperbolicParabolic{typeof(mesh), typeof(equations), typeof(equations_parabolic),
                                         typeof(initial_condition), typeof(_boundary_conditions), typeof(_boundary_conditions_parabolic),
-                                        typeof(source_terms), typeof(solver), typeof(cache), typeof(cache_parabolic)}(
+                                        typeof(source_terms), typeof(solver), typeof(solver_parabolic), typeof(cache), typeof(cache_parabolic)}(
     mesh, equations, equations_parabolic, initial_condition,
     _boundary_conditions, _boundary_conditions_parabolic, source_terms,
-    solver, cache, cache_parabolic)
+    solver, solver_parabolic, cache, cache_parabolic)
 end
 
 
@@ -122,6 +125,7 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s
                                                              equations_parabolic=semi.equations_parabolic,
                                                              initial_condition=semi.initial_condition,
                                                              solver=semi.solver,
+                                                             solver_parabolic=semi.solver_parabolic,
                                                              source_terms=semi.source_terms,
                                                              boundary_conditions=semi.boundary_conditions,
                                                              boundary_conditions_parabolic=semi.boundary_conditions_parabolic
@@ -130,7 +134,7 @@ function remake(semi::SemidiscretizationHyperbolicParabolic; uEltype=real(semi.s
   #       special care if shock-capturing volume integrals are used (because of
   #       the indicators and their own caches...).
   SemidiscretizationHyperbolicParabolic(
-    mesh, equations, equations_parabolic, initial_condition, solver; source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype)
+    mesh, equations, equations_parabolic, initial_condition, solver; solver_parabolic, source_terms, boundary_conditions, boundary_conditions_parabolic, uEltype)
 end
 
 function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic)
@@ -145,6 +149,7 @@ function Base.show(io::IO, semi::SemidiscretizationHyperbolicParabolic)
   print(io, ", ", semi.boundary_conditions_parabolic)
   print(io, ", ", semi.source_terms)
   print(io, ", ", semi.solver)
+  print(io, ", ", semi.solver_parabolic)
   print(io, ", cache(")
   for (idx,key) in enumerate(keys(semi.cache))
     idx > 1 && print(io, " ")
@@ -170,6 +175,7 @@ function Base.show(io::IO, ::MIME"text/plain", semi::SemidiscretizationHyperboli
 
     summary_line(io, "source terms", semi.source_terms)
     summary_line(io, "solver", semi.solver |> typeof |> nameof)
+    summary_line(io, "parabolic solver", semi.solver_parabolic |> typeof |> nameof)
     summary_line(io, "total #DOFs", ndofs(semi))
     summary_footer(io)
   end
@@ -236,7 +242,7 @@ function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t)
 end
 
 function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t)
-  @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, cache, cache_parabolic = semi
+  @unpack mesh, equations_parabolic, initial_condition, boundary_conditions_parabolic, source_terms, solver, solver_parabolic, cache, cache_parabolic = semi
 
   u  = wrap_array(u_ode,  mesh, equations_parabolic, solver, cache_parabolic)
   du = wrap_array(du_ode, mesh, equations_parabolic, solver, cache_parabolic)
@@ -245,7 +251,7 @@ function rhs_parabolic!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabol
   time_start = time_ns()
   @trixi_timeit timer() "parabolic rhs!" rhs_parabolic!(du, u, t, mesh, equations_parabolic, initial_condition,
                                                         boundary_conditions_parabolic, source_terms,
-                                                        solver, cache, cache_parabolic)
+                                                        solver, solver_parabolic, cache, cache_parabolic)
   runtime = time_ns() - time_start
   put!(semi.performance_counter, runtime)
 

diff --git a/src/solvers/dgmulti/dg_parabolic.jl b/src/solvers/dgmulti/dg_parabolic.jl
@@ -1,5 +1,5 @@
 function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsParabolic,
-                                dg::DGMulti, RealT, uEltype)
+                                dg::DGMulti, dg_parabolic, RealT, uEltype)
   nvars = nvariables(equations)
 
   @unpack M, Drst = dg.basis
@@ -19,8 +19,17 @@ function create_cache_parabolic(mesh::DGMultiMesh, equations::AbstractEquationsP
   local_viscous_flux_threaded = [ntuple(_ -> similar(u_transformed, dg.basis.Nq), ndims(mesh)) for _ in 1:Threads.nthreads()]
   local_flux_face_values_threaded = [similar(scalar_flux_face_values[:, 1]) for _ in 1:Threads.nthreads()]
 
+  # precompute 1 / h for penalty terms
+  inv_h = similar(mesh.md.Jf)
+  J = dg.basis.Vf * mesh.md.J # interp to face nodes
+  for e in eachelement(mesh, dg)
+    for i in each_face_node(mesh, dg)
+      inv_h[i, e] = mesh.md.Jf[i, e] / J[i, e]
+    end
+  end
+
   return (; u_transformed, u_grad, viscous_flux,
-            weak_differentiation_matrices,
+            weak_differentiation_matrices, inv_h,
             u_face_values, grad_u_face_values, scalar_flux_face_values,
             local_u_values_threaded, local_viscous_flux_threaded, local_flux_face_values_threaded)
 end
@@ -195,9 +204,25 @@ function calc_viscous_fluxes!(viscous_flux, u, u_grad, mesh::DGMultiMesh,
   end
 end
 
+function calc_viscous_penalty!(scalar_flux_face_values, u_face_values, t, boundary_conditions,
+                               mesh, equations::AbstractEquationsParabolic, dg::DGMulti,
+                               dg_parabolic, cache, cache_parabolic)
+  # compute fluxes at interfaces
+  @unpack scalar_flux_face_values, inv_h = cache_parabolic
+  @unpack mapM, mapP = mesh.md
+  @threaded for face_node_index in each_face_node_global(mesh, dg)
+    idM, idP = mapM[face_node_index], mapP[face_node_index]
+    uM, uP = u_face_values[idM], u_face_values[idP]
+    inv_h_face = inv_h[face_node_index]
+    scalar_flux_face_values[idM] = scalar_flux_face_values[idM] + penalty(uP, uM, inv_h_face, equations, dg_parabolic)
+  end
+  return nothing
+end
+
+
 function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh,
                           equations::AbstractEquationsParabolic,
-                          boundary_conditions, dg::DGMulti, cache, cache_parabolic)
+                          boundary_conditions, dg::DGMulti, dg_parabolic, cache, cache_parabolic)
 
   @unpack weak_differentiation_matrices = cache_parabolic
 
@@ -240,6 +265,10 @@ function calc_divergence!(du, u::StructArray, t, viscous_flux, mesh::DGMultiMesh
   calc_boundary_flux!(scalar_flux_face_values, nothing, t, Divergence(),
                       boundary_conditions, mesh, equations, dg, cache, cache_parabolic)
 
+  calc_viscous_penalty!(scalar_flux_face_values, cache_parabolic.u_face_values, t,
+                        boundary_conditions, mesh, equations, dg, dg_parabolic,
+                        cache, cache_parabolic)
+
   # surface contributions
   apply_to_each_field(mul_by_accum!(dg.basis.LIFT), du, scalar_flux_face_values)
 
@@ -254,7 +283,7 @@ end
 # boundary conditions will be applied to both grad(u) and div(u).
 function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::AbstractEquationsParabolic,
                         initial_condition, boundary_conditions, source_terms,
-                        dg::DGMulti, cache, cache_parabolic)
+                        dg::DGMulti, dg_parabolic, cache, cache_parabolic)
 
   reset_du!(du, dg)
 
@@ -268,7 +297,7 @@ function rhs_parabolic!(du, u, t, mesh::DGMultiMesh, equations_parabolic::Abstra
                        mesh, equations_parabolic, dg, cache, cache_parabolic)
 
   calc_divergence!(du, u_transformed, t, viscous_flux, mesh, equations_parabolic,
-                   boundary_conditions, dg, cache, cache_parabolic)
+                   boundary_conditions, dg, dg_parabolic, cache, cache_parabolic)
 
   return nothing
 

diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl
@@ -7,6 +7,7 @@
 
 include("dg.jl")
 include("dgmulti.jl")
+include("solvers_parabolic.jl")
 
 
 end # @muladd
diff --git a/src/solvers/solvers_parabolic.jl b/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()
diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl
@@ -65,8 +65,8 @@ isdir(outdir) && rm(outdir, recursive=true)
     @test u_flux[2] ≈ u_grad[2]
 
     du = similar(ode.u0)
-    Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic,
-                           boundary_condition_periodic, dg, cache, cache_parabolic)
+    Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic, boundary_condition_periodic,
+                           dg, semi.solver_parabolic, cache, cache_parabolic)
     @test getindex.(du, 1) ≈ 2 * y
   end