Implement BR1 scheme

examples/2d/elixir_advdiff_amr.jl

@@ -0,0 +1,72 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (1.8, 0.6)
+nu = 2.4e-1
+equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+polydeg = 3
+solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs)
+
+coordinates_min = (-4, -4)
+coordinates_max = ( 4,  4)
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_gauss
+semi = SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver_hyperbolic)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=3,
+                                      med_level=4, med_threshold=0.1,
+                                      max_level=5, max_threshold=0.5)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+# stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()

examples/2d/elixir_advdiff_basic.jl

@@ -0,0 +1,68 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advectionvelocity = (0.2, 0.6)
+nu = 1.2e-2
+equations = LinearAdvectionDiffusionEquation2D(advectionvelocity, nu)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+polydeg = 4
+solver_hyperbolic = DGSEM(polydeg, flux_lax_friedrichs)
+
+coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1,  1) # maximum coordinates (max(x), max(y))
+refinement_patches = (
+  (type="box", coordinates_min=(0.0, -1.0), coordinates_max=(1.0, 1.0)),
+)
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                refinement_patches=refinement_patches,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_convergence_test
+semi = SemidiscretizationHyperbolicParabolicBR1(mesh, equations, initial_condition, solver_hyperbolic)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+# stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()

examples/2d/elixir_heat_amr.jl

@@ -0,0 +1,71 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+nu = 4.8e-1
+equations = HeatEquation2D(nu)
+
+# Create DG solver with polynomial degree = 3 and the central flux as surface flux
+solver = DGSEM(3, flux_central)
+
+coordinates_min = (-4, -4) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 4,  4) # maximum coordinates (max(x), max(y))
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=3,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_gauss
+semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first),
+                                      base_level=3,
+                                      med_level=4, med_threshold=0.05,
+                                      max_level=5, max_threshold=0.3)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval=5,
+                           adapt_initial_condition=true,
+                           adapt_initial_condition_only_refine=true)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+# stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, amr_callback)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()

examples/2d/elixir_heat_basic.jl

@@ -0,0 +1,62 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+nu = 1.2e-2
+equations = HeatEquation2D(nu)
+
+# Create DG solver with polynomial degree = 3 and the central flux as surface flux
+solver = DGSEM(3, flux_central)
+
+coordinates_min = (-1, -1) # minimum coordinates (min(x), min(y))
+coordinates_max = ( 1,  1) # maximum coordinates (max(x), max(y))
+
+# Create a uniformely refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=2,
+                n_cells_max=30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+initial_condition = initial_condition_convergence_test
+semi = SemidiscretizationParabolicAuxVars(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 1.0
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100,
+                                     solution_variables=cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+# stepsize_callback = StepsizeCallback(cfl=1.6)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+# callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution)
+
+
+###############################################################################
+# run the simulation
+
+# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+            dt=1e-3, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks, maxiters=1e5);
+
+# Print the timer summary
+summary_callback()

src/Trixi.jl

@@ -59,6 +59,8 @@ include("semidiscretization/semidiscretization_hyperbolic.jl")
 include("callbacks_step/callbacks_step.jl")
 include("callbacks_stage/callbacks_stage.jl")
 include("semidiscretization/semidiscretization_euler_gravity.jl")
+include("semidiscretization/semidiscretization_parabolic_auxvars.jl")
+include("semidiscretization/semidiscretization_hyperbolic_parabolic.jl")
 include("time_integration/time_integration.jl")
 
 # `trixi_include` and special elixirs such as `convergence_test`
@@ -74,7 +76,9 @@ export CompressibleEulerEquations1D, CompressibleEulerEquations2D, CompressibleE
        IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D,
        HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D,
        LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, LinearScalarAdvectionEquation3D,
-       LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D
+       LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
+       HeatEquation2D,
+       LinearAdvectionDiffusionEquation2D
 
 export flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_upwind,
        flux_chandrashekar, flux_chandrashekar_stable, flux_ranocha, flux_derigs_etal, flux_kennedy_gruber, flux_shima_etal
@@ -125,7 +129,9 @@ export DG,
 export nelements, nnodes, nvariables,
        eachelement, eachnode, eachvariable
 
-export SemidiscretizationHyperbolic, semidiscretize, compute_coefficients, integrate
+export SemidiscretizationHyperbolic, SemidiscretizationParabolicAuxVars,
+       SemidiscretizationHyperbolicParabolic, SemidiscretizationHyperbolicParabolicBR1,
+       semidiscretize, compute_coefficients, integrate
 
 export SemidiscretizationEulerGravity, ParametersEulerGravity,
        timestep_gravity_erk52_3Sstar!, timestep_gravity_carpenter_kennedy_erk54_2N!

src/equations/equations.jl

@@ -55,10 +55,12 @@ default_analysis_integrals(::AbstractEquations)  = (entropy_timederivative,)
 
 """
     flux_central(u_ll, u_rr, orientation, equations::AbstractEquations)
+    flux_central(u_ll, u_rr, gradients_ll, gradients_rr, orientation, equations::AbstractEquations)
 
 The classical central numerical flux `f((u_ll) + f(u_rr)) / 2`. When this flux is
 used as volume flux, the discretization is equivalent to the classical weak form
-DG method (except floating point errors).
+DG method (except floating point errors). The second version call the parabolic flux, which requires
+the gradients as additional arguments.
 """
 @inline function flux_central(u_ll, u_rr, orientation, equations::AbstractEquations)
   # Calculate regular 1D fluxes
@@ -68,6 +70,14 @@ DG method (except floating point errors).
   # Average regular fluxes
   return 0.5 * (f_ll + f_rr)
 end
+@inline function flux_central(u_ll, u_rr, gradients_ll, gradients_rr, orientation, equations::AbstractEquations)
+  # Calculate regular 1D fluxes
+  f_ll = calcflux(u_ll, gradients_ll, orientation, equations)
+  f_rr = calcflux(u_rr, gradients_rr, orientation, equations)
+
+  # Average regular fluxes
+  return 0.5 * (f_ll + f_rr)
+end
 
 
 @inline cons2cons(u, ::AbstractEquations) = u
@@ -114,3 +124,15 @@ include("hyperbolic_diffusion_3d.jl")
 abstract type AbstractLatticeBoltzmannEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
 include("lattice_boltzmann_2d.jl")
 include("lattice_boltzmann_3d.jl")
+
+# Gradient equations
+abstract type AbstractGradientEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+include("gradient_equations_2d.jl")
+
+# Heat equation
+abstract type AbstractHeatEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+include("heat_equation_2d.jl")
+
+# Linear scalar advection-diffusion equation
+abstract type AbstractLinearAdvectionDiffusionEquation{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end
+include("linear_advection_diffusion_2d.jl")

src/equations/gradient_equations_2d.jl

@@ -0,0 +1,46 @@
+
+@doc raw"""
+    GradientEquations2D
+
+The gradient equations
+```math
+q^d - \partial_d u = 0
+```
+in direction `d` in two space dimensions as required for, e.g., the Bassi & Rebay 1 (BR1) or the
+local discontinuous Galerkin (LDG) schemes.
+"""
+struct GradientEquations2D{RealT<:Real, NVARS} <: AbstractGradientEquations{2, NVARS}
+  orientation::Int
+end
+
+GradientEquations2D(::Type{RealT}, nvars, orientation) where RealT = GradientEquations2D{RealT, nvars}(orientation)
+GradientEquations2D(nvars, orientation) = GradientEquations2D(Float64, nvars, orientation)
+
+
+get_name(::GradientEquations2D) = "GradientEquations2D"
+varnames(::typeof(cons2cons), equations::GradientEquations2D) = SVector(ntuple(v -> "gradient_"*string(v), nvariables(equations)))
+varnames(::typeof(cons2prim), equations::GradientEquations2D) = varnames(cons2cons, equations)
+
+# Set initial conditions at physical location `x` for time `t`
+"""
+    initial_condition_constant(x, t, equations::GradientEquations2D)
+
+A constant initial condition to test free-stream preservation.
+"""
+function initial_condition_constant(x, t, equations::GradientEquations2D)
+  return SVector(ntuple(v -> zero(eltype(x)), nvariables(equations)))
+end
+
+
+# Pre-defined source terms should be implemented as
+# function source_terms_WHATEVER(u, x, t, equations::GradientEquations2D)
+
+
+# Calculate 1D flux in for a single point
+@inline function calcflux(u, orientation, equations::GradientEquations2D)
+  if orientation == equations.orientation
+    return -u
+  else
+    return SVector(ntuple(v -> zero(eltype(u)), nvariables(equations)))
+  end
+end