Parabolic terms for DGMulti (#1136)

* minor formatting * adding multi-term semidiscretization * fixing bug and adding test * changing name * changing name, importing SplitODEFunction * fixing test by adding norm * trying another test fix * adding norm back * adjusting tol * try to fix test again * syntax error * removing norm frm test * adding cache to parabolic terms (enable reuse of inviscid cache) * generalize dgmulti iterators * fix typo * generalizing reset_du! for dgmulti * WIP gradient and divergence computations * fix name * more generalization * divergence + remove repeated code * adding cache * remove test * update diffusion test * update diffusion rhs * moving routines around * fix multiple includes * removing time dependence for now from viscous rhs * export ScalarDiffusion2D * fix test * working periodic advection RHS * cleanup and comments * making periodic advection-diffusion example * trim whitespace * adding t back in to parabolic rhs add t back in * refactoring boundary condition code * update elixir with BCs * change default DGMulti interface flux behavior same as in #1139 * cleanup * renaming variables * BC prepping * cleanup of elixir * more name fixing * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Hendrik Ranocha <[email protected]> * adding Neumann BCs * switching to weak form weak form more weak form weak form * fixing cache variable names * simplified parabolic boundary flux treatement * updated elixir with different types of BCs * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Update src/equations/equations.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * rhs! -> rhs_parabolic! * move Gradient/Divergence types * add "no boundary condition" BC * rename ScalarDiffusion -> LaplaceDiffusion r r * rhs -> rhs_parabolic! again * comments * moving tests around * scalarDiffusion -> LaplaceDiffusion * tuple-based constructor * updating examples * add ScalarPlotData2D docstring (#1145) Co-authored-by: Jesse Chan <[email protected]> * Apply suggestions from code review Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]> * renaming for consistency parabolic_equations -> equations_parabolic no_boundary_condition -> boundary_condition_do_nothing * move Neumann BCs into equations.jl * renaming, e.g., ParabolicEquations -> EquationsParabolic, etc * rename test module * generalize LaplaceDiffusion2D to any number of variables * removing saveat times and standardizing `tol` names * parabolic_cache -> cache_parabolic cache * renaming `create_cache` -> `create_cache_parabolic` avoids conflicts when using dg::DGMultiFluxDiff * Update examples/dgmulti_2d/elixir_advection_diffusion.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Update src/equations/equations.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * parabolic_boundary_conditions -> boundary_conditions_parabolic similarly for hyperbolic d * Update src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl Co-authored-by: Michael Schlottke-Lakemper <[email protected]> * add equations as field to LaplaceDiffusion2D * moving Neumann BCs into equations.jl and adding default BC behavior * running parabolic tests, adding integration test * adding parabolic tests to ci * missing comma * Update src/equations/laplace_diffusion_2d.jl Co-authored-by: Hendrik Ranocha <[email protected]> * fix test * fix test typo * fix tests (for real this time) * improve test coverage Co-authored-by: Jesse Chan <[email protected]> Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
trixi-framework · May 27, 2022 · c60c025 · c60c025
1 parent b6d801c
commit c60c025
Show file tree

Hide file tree

Showing 14 changed files with 822 additions and 21 deletions.
diff --git a/.github/workflows/ci.yml b/.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

diff --git a/examples/dgmulti_2d/elixir_advection_diffusion.jl b/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(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
+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
diff --git a/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl b/examples/dgmulti_2d/elixir_advection_diffusion_periodic.jl
@@ -0,0 +1,36 @@
+
+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, 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, 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
+
+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
diff --git a/src/Trixi.jl b/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,

diff --git a/src/equations/equations.jl b/src/equations/equations.jl
@@ -175,6 +175,45 @@ 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 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
 """
@@ -329,5 +368,10 @@ 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
+
+# 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/laplace_diffusion_2d.jl b/src/equations/laplace_diffusion_2d.jl
@@ -0,0 +1,32 @@
+@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`.
+"""
+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 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::BoundaryConditionNeumann)(flux_inner, normal::AbstractVector,
+                                                                x, t, operator_type::Divergence,
+                                                                equations::LaplaceDiffusion2D)
+  return boundary_condition.boundary_normal_flux_function(x, t, equations)
+end