Incorporate parabolic terms into main

.github/workflows/ci.yml

@@ -69,6 +69,7 @@ jobs:
           - p4est_part1
           - p4est_part2
           - unstructured_dgmulti
+          - parabolic
           - paper_self_gravitating_gas_dynamics
           - misc_part1
           - misc_part2

NEWS.md

@@ -9,6 +9,22 @@ for human readability.
 
 #### Added
 
+- Experimental support for 2D parabolic diffusion terms has been added.
+  * `LaplaceDiffusion2D` and `CompressibleNavierStokesDiffusion2D` can be used to add
+  diffusion to systems. `LaplaceDiffusion2D` can be used to add scalar diffusion to each
+  equation of a system, while `CompressibleNavierStokesDiffusion2D` can be used to add
+  Navier-Stokes diffusion to `CompressibleEulerEquations2D`.
+  * Parabolic boundary conditions can be imposed as well. For `LaplaceDiffusion2D`, both
+  `Dirichlet` and `Neumann` conditions are supported. For `CompressibleNavierStokesDiffusion2D`,
+  viscous no-slip velocity boundary conditions are supported, along with adiabatic and isothermal
+  temperature boundary conditions. See the boundary condition container
+  `BoundaryConditionNavierStokesWall` and boundary condition types `NoSlip`, `Adiabatic`, and
+  `Isothermal` for more information.
+  * `CompressibleNavierStokesDiffusion2D` can utilize both primitive variables (which are not
+  guaranteed to provably dissipate entropy) and entropy variables (which provably dissipate
+  entropy at the semi-discrete level).
+  * Please check the `examples` directory for further information about the supported setups.
+    Further documentation will be added later.
 - Numerical fluxes `flux_shima_etal_turbo` and `flux_ranocha_turbo` that are
   equivalent to their non-`_turbo` counterparts but may enable specialized
   methods making use of SIMD instructions to increase runtime efficiency

docs/literate/src/files/adding_new_parabolic_terms.jl

@@ -0,0 +1,160 @@
+#src # Adding new parabolic terms.
+
+# This demo illustrates the steps involved in adding new parabolic terms for the scalar
+# advection equation. In particular, we will add an anisotropic diffusion. We begin by
+# defining the hyperbolic (advection) part of the advection-diffusion equation.
+
+using OrdinaryDiffEq
+using Trixi
+
+
+advection_velocity = (1.0, 1.0)
+equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity);
+
+# ## Define a new parabolic equation type
+#
+# Next, we define a 2D parabolic diffusion term type. This is similar to [`LaplaceDiffusion2D`](@ref)
+# except that the `diffusivity` field refers to a spatially constant diffusivity matrix now. Note that
+# `ConstantAnisotropicDiffusion2D` has a field for `equations_hyperbolic`. It is useful to have
+# information about the hyperbolic system available to the parabolic part so that we can reuse
+# functions defined for hyperbolic equations (such as `varnames`).
+
+struct ConstantAnisotropicDiffusion2D{E, T} <: Trixi.AbstractEquationsParabolic{2, 1}
+  diffusivity::T
+  equations_hyperbolic::E
+end
+
+varnames(variable_mapping, equations_parabolic::ConstantAnisotropicDiffusion2D) =
+  varnames(variable_mapping, equations_parabolic.equations_hyperbolic)
+
+# Next, we define the viscous flux function. We assume that the mixed hyperbolic-parabolic system
+# is of the form
+# ```math
+# \partial_t u(t,x) + \partial_x (f_1(u) - g_1(u, \nabla u))
+#                   + \partial_y (f_2(u) - g_2(u, \nabla u)) = 0
+# ```
+# where ``f_1(u)``, ``f_2(u)`` are the hyperbolic fluxes and ``g_1(u, \nabla u)``, ``g_2(u, \nabla u)`` denote
+# the viscous fluxes. For anisotropic diffusion, the viscous fluxes are the first and second components
+# of the matrix-vector product involving `diffusivity` and the gradient vector.
+#
+# Here, we specialize the flux to our new parabolic equation type `ConstantAnisotropicDiffusion2D`.
+
+function Trixi.flux(u, gradients, orientation::Integer, equations_parabolic::ConstantAnisotropicDiffusion2D)
+  @unpack diffusivity = equations_parabolic
+  dudx, dudy = gradients
+  if orientation == 1
+    return SVector(diffusivity[1, 1] * dudx + diffusivity[1, 2] * dudy)
+  else # if orientation == 2
+    return SVector(diffusivity[2, 1] * dudx + diffusivity[2, 2] * dudy)
+  end
+end
+
+# ## Defining boundary conditions
+
+# Trixi.jl's implementation of parabolic terms discretizes both the gradient and divergence
+# using weak formulation. In other words, we discretize the system
+# ```math
+# \begin{aligned}
+# \bm{q} &= \nabla u \\
+# \bm{\sigma} &= \begin{pmatrix} g_1(u, \bm{q}) \\ g_2(u, \bm{q}) \end{pmatrix} \\
+# \text{viscous contribution } &= \nabla \cdot \bm{\sigma}
+# \end{aligned}
+# ```
+#
+# Boundary data must be specified for all spatial derivatives, e.g., for both the gradient
+# equation ``\bm{q} = \nabla u`` and the divergence of the viscous flux
+# ``\nabla \cdot \bm{\sigma}``. We account for this by introducing internal `Gradient`
+# and `Divergence` types which are used to dispatch on each type of boundary condition.
+#
+# As an example, let us introduce a Dirichlet boundary condition with constant boundary data.
+
+struct BoundaryConditionConstantDirichlet{T <: Real}
+  boundary_value::T
+end
+
+# This boundary condition contains only the field `boundary_value`, which we assume to be some
+# real-valued constant which we will impose as the Dirichlet data on the boundary.
+#
+# Boundary conditions have generally been defined as "callable structs" (also known as "functors").
+# For each boundary condition, we need to specify the appropriate boundary data to return for both
+# the `Gradient` and `Divergence`. Since the gradient is operating on the solution `u`, the boundary
+# data should be the value of `u`, and we can directly impose Dirichlet data.
+
+@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector,
+                                                                          x, t, operator_type::Trixi.Gradient,
+                                                                          equations_parabolic::ConstantAnisotropicDiffusion2D)
+  return boundary_condition.boundary_value
+end
+
+# While the gradient acts on the solution `u`, the divergence acts on the viscous flux ``\bm{\sigma}``.
+# Thus, we have to supply boundary data for the `Divergence` operator that corresponds to ``\bm{\sigma}``.
+# However, we've already imposed boundary data on `u` for a Dirichlet boundary condition, and imposing
+# boundary data for ``\bm{\sigma}`` might overconstrain our problem.
+#
+# Thus, for the `Divergence` boundary data under a Dirichlet boundary condition, we simply return
+# `flux_inner`, which is boundary data for ``\bm{\sigma}`` computed using the "inner" or interior solution.
+# This way, we supply boundary data for the divergence operation without imposing any additional conditions.
+
+@inline function (boundary_condition::BoundaryConditionConstantDirichlet)(flux_inner, u_inner, normal::AbstractVector,
+                                                                          x, t, operator_type::Trixi.Divergence,
+                                                                          equations_parabolic::ConstantAnisotropicDiffusion2D)
+  return flux_inner
+end
+
+# ### A note on the choice of gradient variables
+#
+# It is often simpler to transform the solution variables (and solution gradients) to another set of
+# variables prior to computing the viscous fluxes (see [`CompressibleNavierStokesDiffusion2D`](@ref)
+# for an example of this). If this is done, then the boundary condition for the `Gradient` operator
+# should be modified accordingly as well.
+#
+# ## Putting things together
+#
+# Finally, we can instantiate our new parabolic equation type, define boundary conditions,
+# and run a simulation. The specific anisotropic diffusion matrix we use produces more
+# dissipation in the direction ``(1, -1)`` as an isotropic diffusion.
+#
+# For boundary conditions, we impose that ``u=1`` on the left wall, ``u=2`` on the bottom
+# wall, and ``u = 0`` on the outflow walls. The initial condition is taken to be ``u = 0``.
+# Note that we use `BoundaryConditionConstantDirichlet` only for the parabolic boundary
+# conditions, since we have not defined its behavior for the hyperbolic part.
+
+using Trixi: SMatrix
+diffusivity = 5.0e-2 * SMatrix{2, 2}([2 -1; -1 2])
+equations_parabolic = ConstantAnisotropicDiffusion2D(diffusivity, equations_hyperbolic);
+
+boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1.0)),
+                                    y_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(2.0)),
+                                    y_pos = boundary_condition_do_nothing,
+                                    x_pos = boundary_condition_do_nothing)
+
+boundary_conditions_parabolic = (; x_neg = BoundaryConditionConstantDirichlet(1.0),
+                                   y_neg = BoundaryConditionConstantDirichlet(2.0),
+                                   y_pos = BoundaryConditionConstantDirichlet(0.0),
+                                   x_pos = BoundaryConditionConstantDirichlet(0.0));
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                periodicity=false, n_cells_max=30_000) # set maximum capacity of tree data structure
+
+initial_condition = (x, t, equations) -> SVector(0.0)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations_hyperbolic, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions_hyperbolic,
+                                                                  boundary_conditions_parabolic))
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+callbacks = CallbackSet(SummaryCallback())
+time_int_tol = 1.0e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            save_everystep=false, callback=callbacks);
+
+using Plots
+plot(sol)
+

docs/literate/src/files/parabolic_terms.jl

@@ -0,0 +1,88 @@
+#src # Parabolic terms (advection-diffusion).
+
+# Experimental support for parabolic diffusion terms is available in Trixi.jl.
+# This demo illustrates parabolic terms for the advection-diffusion equation.
+
+using OrdinaryDiffEq
+using Trixi
+
+# ## Splitting a system into hyperbolic and parabolic parts.
+
+# For a mixed hyperbolic-parabolic system, we represent the hyperbolic and parabolic
+# parts of the system  separately. We first define the hyperbolic (advection) part of
+# the advection-diffusion equation.
+
+advection_velocity = (1.5, 1.0)
+equations_hyperbolic = LinearScalarAdvectionEquation2D(advection_velocity);
+
+# Next, we define the parabolic diffusion term. The constructor requires knowledge of
+# `equations_hyperbolic` to be passed in because the [`LaplaceDiffusion2D`](@ref) applies
+# diffusion to every variable of the hyperbolic system.
+
+diffusivity = 5.0e-2
+equations_parabolic = LaplaceDiffusion2D(diffusivity, equations_hyperbolic);
+
+# ## Boundary conditions
+
+# As with the equations, we define boundary conditions separately for the hyperbolic and
+# parabolic part of the system. For this example, we impose inflow BCs for the hyperbolic
+# system (no condition is imposed on the outflow), and we impose Dirichlet boundary conditions
+# for the parabolic equations. Both `BoundaryConditionDirichlet` and `BoundaryConditionNeumann`
+# are defined for `LaplaceDiffusion2D`.
+#
+# The hyperbolic and parabolic boundary conditions are assumed to be consistent with each other.
+
+boundary_condition_zero_dirichlet = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))
+
+boundary_conditions_hyperbolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])),
+                                    y_neg = boundary_condition_zero_dirichlet,
+                                    y_pos = boundary_condition_do_nothing,
+                                    x_pos = boundary_condition_do_nothing)
+
+boundary_conditions_parabolic = (; x_neg = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.5 * x[2])),
+                                   y_neg = boundary_condition_zero_dirichlet,
+                                   y_pos = boundary_condition_zero_dirichlet,
+                                   x_pos = boundary_condition_zero_dirichlet);
+
+# ## Defining the solver and mesh
+
+# The process of creating the DG solver and mesh is the same as for a purely
+# hyperbolic system of equations.
+
+solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)
+coordinates_min = (-1.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1.0,  1.0) # maximum coordinates (max(x), max(y))
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                periodicity=false, n_cells_max=30_000) # set maximum capacity of tree data structure
+
+initial_condition = (x, t, equations) -> SVector(0.0);
+
+# ## Semidiscretizing and solving
+
+# To semidiscretize a hyperbolic-parabolic system, we create a [`SemidiscretizationHyperbolicParabolic`](@ref).
+# This differs from a [`SemidiscretizationHyperbolic`](@ref) in that we pass in a `Tuple` containing both the
+# hyperbolic and parabolic equation, as well as a `Tuple` containing the hyperbolic and parabolic
+# boundary conditions.
+
+semi = SemidiscretizationHyperbolicParabolic(mesh,
+                                             (equations_hyperbolic, equations_parabolic),
+                                             initial_condition, solver;
+                                             boundary_conditions=(boundary_conditions_hyperbolic,
+                                                                  boundary_conditions_parabolic))
+
+# The rest of the code is identical to the hyperbolic case. We create a system of ODEs through
+# `semidiscretize`, defining callbacks, and then passing the system to OrdinaryDiffEq.jl.
+
+tspan = (0.0, 1.5)
+ode = semidiscretize(semi, tspan)
+callbacks = CallbackSet(SummaryCallback())
+time_int_tol = 1.0e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            save_everystep=false, callback=callbacks);
+
+# We can now visualize the solution, which develops a boundary layer at the outflow boundaries.
+
+using Plots
+plot(sol)
+

docs/make.jl

@@ -39,6 +39,8 @@ files = [
     # Topic: equations
     "Adding a new scalar conservation law" => "adding_new_scalar_equations.jl",
     "Adding a non-conservative equation" => "adding_nonconservative_equation.jl",
+    "Parabolic terms" => "parabolic_terms.jl",
+    "Adding new parabolic terms" => "adding_new_parabolic_terms.jl",
     # Topic: meshes
     "Adaptive mesh refinement" => "adaptive_mesh_refinement.jl",
     "Structured mesh with curvilinear mapping" => "structured_mesh_mapping.jl",

examples/dgmulti_2d/elixir_advection_diffusion.jl

@@ -0,0 +1,56 @@
+using Trixi, OrdinaryDiffEq
+
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
+             volume_integral = VolumeIntegralWeakForm())
+
+equations = LinearScalarAdvectionEquation2D(1.5, 1.0)
+equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations)
+
+initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0)
+initial_condition = initial_condition_zero
+
+# tag different boundary segments
+left(x, tol=50*eps()) = abs(x[1] + 1) < tol
+right(x, tol=50*eps()) = abs(x[1] - 1) < tol
+bottom(x, tol=50*eps()) = abs(x[2] + 1) < tol
+top(x, tol=50*eps()) = abs(x[2] - 1) < tol
+is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom)
+mesh = DGMultiMesh(dg, cells_per_dimension=(16, 16); is_on_boundary)
+
+# BC types
+boundary_condition_left = BoundaryConditionDirichlet((x, t, equations) -> SVector(1 + 0.1 * x[2]))
+boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations) -> SVector(0.0))
+boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations) -> SVector(0.0))
+
+# define inviscid boundary conditions
+boundary_conditions = (; :left => boundary_condition_left,
+                         :bottom => boundary_condition_zero,
+                         :top => boundary_condition_do_nothing,
+                         :right => boundary_condition_do_nothing)
+
+# define viscous boundary conditions
+boundary_conditions_parabolic = (; :left => boundary_condition_left,
+                                   :bottom => boundary_condition_zero,
+                                   :top => boundary_condition_zero,
+                                   :right => boundary_condition_neumann_zero)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic), initial_condition, dg;
+                                             boundary_conditions=(boundary_conditions, boundary_conditions_parabolic))
+
+tspan = (0.0, 1.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
+callbacks = CallbackSet(summary_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol,
+            save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary