Parabolic terms for DGMulti

examples/dgmulti_2d/elixir_advection_diffusion.jl

@@ -0,0 +1,55 @@
+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(5e-2)
+
+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
+parabolic_boundary_conditions = (; :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, parabolic_boundary_conditions))
+
+tspan = (0.0, 1.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval=10)
+callbacks = CallbackSet(summary_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-6
+saveat = LinRange(tspan..., 10)
+sol = solve(ode, RDPK3SpFSAL49(), abstol=time_int_tol, reltol=time_int_tol, save_everystep=false,
+            saveat=saveat, callback=callbacks)
+summary_callback() # print the timer summary

examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl

@@ -0,0 +1,34 @@
+
+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.0, 1.0)
+equations_parabolic = LaplaceDiffusion2D(1e-2)
+
+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, equations_parabolic),
+                                             initial_condition, dg)
+
+tspan = (0.0, 2.0)
+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, analysis_callback)
+
+###############################################################################
+# run the simulation
+
+tol = 1e-6
+sol = solve(ode, RDPK3SpFSAL49(), abstol=tol, reltol=tol, save_everystep=false, callback=callbacks)
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -26,7 +26,8 @@ using SparseArrays: AbstractSparseMatrix, AbstractSparseMatrixCSC, sparse, dropt
 using Reexport: @reexport
 
 using SciMLBase: CallbackSet, DiscreteCallback,
-                 ODEProblem, ODESolution, ODEFunction
+                 ODEProblem, ODESolution, ODEFunction,
+                 SplitODEProblem
 import SciMLBase: get_du, get_tmp_cache, u_modified!,
                   AbstractODEIntegrator, init, step!, check_error,
                   get_proposed_dt, set_proposed_dt!,
@@ -106,6 +107,7 @@ include("meshes/meshes.jl")
 include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization.jl")
 include("semidiscretization/semidiscretization_hyperbolic.jl")
+include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("semidiscretization/semidiscretization_euler_acoustics.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
@@ -128,6 +130,7 @@ export AcousticPerturbationEquations2D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D,
        LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D,
        InviscidBurgersEquation1D,
+       LaplaceDiffusion2D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D
 
@@ -151,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
@@ -185,6 +190,8 @@ export nelements, nnodes, nvariables,
 
 export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate
 
+export SemidiscretizationHyperbolicParabolic
+
 export SemidiscretizationEulerAcoustics
 
 export SemidiscretizationEulerGravity, ParametersEulerGravity,

src/equations/equations.jl

@@ -175,6 +175,25 @@ end
   return flux
 end
 
+"""
+    BoundaryConditionNeumann(boundary_value_function)
+
+Similar to `BoundaryConditionDirichlet`, but creates a Neumann boundary condition for parabolic
+equations that uses the function  `boundary_value_function` to specify the values 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_value_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_value_function::B
+end
+
+# Imposing no boundary condition just evaluates the flux at the inner state.
+@inline function boundary_condition_do_nothing(flux, other_args...)
+  return flux
+end
 
 # set sensible default values that may be overwritten by specific equations
 """
@@ -329,5 +348,14 @@ include("lattice_boltzmann_3d.jl")
 abstract type AbstractAcousticPerturbationEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
 include("acoustic_perturbation_2d.jl")
 
+abstract type AbstractEquationsParabolic{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+
+# operator types used for dispatch on parabolic boundary fluxes
+struct Gradient end
+struct Divergence 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

src/equations/laplace_diffusion_2d.jl

@@ -0,0 +1,53 @@
+@doc raw"""
+    LaplaceDiffusion2D(diffusivity, equations)
+
+`LaplaceDiffusion2D` represents a scalar diffusion term ``\nabla \cdot (\kappa\nabla u))``
+with diffusivity ``\kappa`` applied to each solution component defined by `equations`.
+
+If `equations` is not specified, a scalar `LaplaceDiffusion2D` is returned.
+"""
+struct LaplaceDiffusion2D{E, N, T} <: AbstractLaplaceDiffusionEquations{2, N}
+  diffusivity::T
+  equations::E
+end
+
+function LaplaceDiffusion2D(diffusivity, equations=nothing)
+  if equations isa AbstractEquations
+    nvars = nvariables(equations)
+  else # assume scalar diffusion if no equations are specified
+    nvars = 1
+  end
+  return LaplaceDiffusion2D{typeof(equations), nvars, typeof(diffusivity)}(diffusivity, equations)
+end
+
+# no orientation specified since the flux is vector-valued
+function flux(u, grad_u, equations::LaplaceDiffusion2D)
+  dudx, dudy = grad_u
+  return SVector(equations.diffusivity * dudx, equations.diffusivity * dudy)
+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,
+                                                                  equations::LaplaceDiffusion2D)
+  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
+
+
+@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)
+end