Parabolic terms for DGMulti

.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

examples/dgmulti_2d/elixir_advection_diffusion.jl

@@ -1,43 +1,42 @@
-
 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)
-parabolic_equations = LaplaceDiffusion2D(5e-2)
+equations_parabolic = LaplaceDiffusion2D(5e-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
+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 + .1 * x[2]))
+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 => no_boundary_condition,
-                         :right => no_boundary_condition)
+                         :top => boundary_condition_do_nothing,
+                         :right => boundary_condition_do_nothing)
 
 # define viscous boundary conditions
-parabolic_boundary_conditions = (; :left => boundary_condition_left,
+boundary_conditions_parabolic = (; :left => boundary_condition_left,
                                    :bottom => boundary_condition_zero,
                                    :top => boundary_condition_zero,
                                    :right => boundary_condition_neumann_zero)
 
-semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, parabolic_equations), initial_condition, dg;
-                                             boundary_conditions=(boundary_conditions, parabolic_boundary_conditions))
+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)
@@ -49,8 +48,7 @@ callbacks = CallbackSet(summary_callback, alive_callback)
 ###############################################################################
 # run the simulation
 
-tol = 1e-6
-tsave = LinRange(tspan..., 10)
-sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false,
-            saveat=tsave, callback=callbacks)
+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

examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl

@@ -6,15 +6,15 @@ dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(
              volume_integral = VolumeIntegralWeakForm())
 
 equations = LinearScalarAdvectionEquation2D(1.0, 1.0)
-parabolic_equations = LaplaceDiffusion2D(1e-2)
+equations_parabolic = LaplaceDiffusion2D(1e-2, equations)
 
 function initial_condition_sharp_gaussian(x, t, equations::LinearScalarAdvectionEquation2D)
   return SVector(exp(-100 * (x[1]^2 + x[2]^2)))
 end
 initial_condition = initial_condition_sharp_gaussian
 
 mesh = DGMultiMesh(dg, cells_per_dimension = (16, 16), periodicity=true)
-semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, parabolic_equations),
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                              initial_condition, dg)
 
 tspan = (0.0, 2.0)
@@ -29,6 +29,8 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 ###############################################################################
 # run the simulation
 
-tol = 1e-6
-sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, callback=callbacks)
+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

src/Trixi.jl

@@ -130,13 +130,10 @@ export AcousticPerturbationEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D,
        LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D,
        InviscidBurgersEquation1D,
+       LaplaceDiffusion2D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D
 
-export no_boundary_condition
-
-export LaplaceDiffusion2D, BoundaryConditionNeumann
-
 export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov,
        flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner,
        flux_nonconservative_powell,
@@ -157,8 +154,10 @@ export initial_condition_constant,
        initial_condition_density_wave,
        initial_condition_weak_blast_wave
 
-export boundary_condition_periodic,
+export boundary_condition_do_nothing,
+       boundary_condition_periodic,
        BoundaryConditionDirichlet,
+       BoundaryConditionNeumann,
        boundary_condition_noslip_wall,
        boundary_condition_slip_wall,
        boundary_condition_wall

src/equations/equations.jl

@@ -175,9 +175,44 @@ end
   return flux
 end
 
+# operator types used for dispatch on parabolic boundary fluxes
+struct Gradient end
+struct Divergence end
+
+# Dirichlet-type boundary condition for use with parabolic equations.
+@inline function (boundary_condition::BoundaryConditionDirichlet)(u_inner, normal::AbstractVector,
+                                                                  x, t, operator_type::Divergence,
+                                                                  equations)
+  return u_inner
+end
+
+
+"""
+    BoundaryConditionNeumann(boundary_normal_flux_function)
+
+Similar to `BoundaryConditionDirichlet`, but creates a Neumann boundary condition for parabolic
+equations that uses the function `boundary_normal_flux_function` to specify the values of the normal
+flux at the boundary.
+The passed boundary value function will be called with the same arguments as an initial condition function is called, i.e., as
+```julia
+boundary_normal_flux_function(x, t, equations)
+```
+where `x` specifies the coordinates, `t` is the current time, and `equation` is the corresponding system of equations.
+"""
+struct BoundaryConditionNeumann{B}
+  boundary_normal_flux_function::B
+end
+
+# specify default behavior for a Neumann BC is to evaluate the flux at the inner state
+@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector,
+                                                                x, t, operator_type::Gradient,
+                                                                equations)
+  return flux_inner
+end
+
 # Imposing no boundary condition just evaluates the flux at the inner state.
-@inline function no_boundary_condition(flux, other_args...)
-  return flux
+@inline function boundary_condition_do_nothing(inner_flux_or_state, other_args...)
+  return inner_flux_or_state
 end
 
 # set sensible default values that may be overwritten by specific equations
@@ -333,14 +368,10 @@ include("lattice_boltzmann_3d.jl")
 abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
 include("acoustic_perturbation_2d.jl")
 
-abstract type AbstractParabolicEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
-
-# operator types used for dispatch on parabolic boundary fluxes
-struct Gradient end
-struct Divergence end
+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} <: AbstractParabolicEquations{NDIMS, NVARS} end
+abstract type AbstractLaplaceDiffusionEquations{NDIMS, NVARS} <: AbstractEquationsParabolic{NDIMS, NVARS} end
 include("laplace_diffusion_2d.jl")
 
 end # @muladd

src/equations/laplace_diffusion_2d.jl

@@ -1,17 +1,23 @@
 @doc raw"""
-    LaplaceDiffusion2D{T} <: AbstractParabolicEquations{2, 1}
+    LaplaceDiffusion2D(diffusivity, equations)
 
 `LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))``
-applied to each solution component.
+with diffusivity ``\kappa`` applied to each solution component defined by `equations`.
+
+If `equations` is not specified, a scalar `LaplaceDiffusion2D` is returned.
 """
-struct LaplaceDiffusion2D{T} <: AbstractLaplaceDiffusionEquations{2, 1}
+struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N}
   diffusivity::T
+  equations::E
 end
 
+LaplaceDiffusion2D(diffusivity, equations) =
+  LaplaceDiffusion2D{typeof(equations), nvariables(equations), typeof(diffusivity)}(diffusivity, equations)
+
 # no orientation specified since the flux is vector-valued
 function flux(u, grad_u, equations::LaplaceDiffusion2D)
   dudx, dudy = grad_u
-  return equations.diffusivity * dudx, equations.diffusivity * dudy
+  return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy)
 end
 
 # Dirichlet-type boundary condition for use with a parabolic solver in weak form
@@ -21,26 +27,8 @@ end
   return boundary_condition.boundary_value_function(x, t, equations)
 end
 
-@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, normal::AbstractVector,
-                                                                  x, t, operator_type::Divergence,
-                                                                  equations::LaplaceDiffusion2D)
-  return flux_inner
-end
-
-# `boundary_value_function` should have signature `boundary_value_function(x, t, equations)`
-# For Neumann BCs, this corresponds to kappa * (∇u ⋅ n)
-struct BoundaryConditionNeumann{B}
-  boundary_value_function::B
-end
-
-@inline function (boundary_condition::BoundaryConditionNeumann)(u_inner, normal::AbstractVector,
-                                                                x, t, operator_type::Gradient,
-                                                                equations::LaplaceDiffusion2D)
-  return u_inner
-end
-
 @inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector,
                                                                 x, t, operator_type::Divergence,
                                                                 equations::LaplaceDiffusion2D)
-  return boundary_condition.boundary_value_function(x, t, equations)
+  return boundary_condition.boundary_normal_flux_function(x, t, equations)
 end