Add Laplace diffusion in terms of entropy variables

examples/dgmulti_2d/elixir_euler_laplace_diffusion.jl

@@ -0,0 +1,38 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+surface_flux = FluxLaxFriedrichs()
+volume_flux = flux_ranocha
+dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(surface_flux),
+             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+equations = CompressibleEulerEquations2D(1.4)
+equations_parabolic = LaplaceDiffusionEntropyVariables2D(0.001, equations)
+
+initial_condition = initial_condition_weak_blast_wave
+
+cells_per_dimension = (16, 16)
+mesh = DGMultiMesh(dg, cells_per_dimension,
+                   coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0),
+                   periodicity = true)
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, dg)
+
+tspan = (0.0, 1.1)
+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,
+                        analysis_callback,
+                        alive_callback)
+
+###############################################################################
+# run the simulation
+
+alg = RDPK3SpFSAL35()
+sol = solve(ode, alg; abstol = 1.0e-6, reltol = 1.0e-6,
+            ode_default_options()..., callback = callbacks);

examples/tree_1d_dgsem/elixir_euler_laplace_diffusion.jl

@@ -0,0 +1,54 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquations1D(1.4)
+equations_parabolic = LaplaceDiffusionEntropyVariables1D(0.01, equations)
+
+initial_condition = initial_condition_weak_blast_wave
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0,)
+coordinates_max = (2.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 5,
+                n_cells_max = 10_000)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

examples/tree_3d_dgsem/elixir_euler_laplace_diffusion.jl

@@ -0,0 +1,55 @@
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+equations_parabolic = LaplaceDiffusionEntropyVariables3D(0.01, equations)
+
+initial_condition = initial_condition_weak_blast_wave
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0, -2.0, -2.0)
+coordinates_max = (2.0, 2.0, 2.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 100_000)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.3)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

src/Trixi.jl

@@ -174,6 +174,8 @@ export AcousticPerturbationEquations2D,
        PassiveTracerEquations
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
+       LaplaceDiffusionEntropyVariables1D, LaplaceDiffusionEntropyVariables2D,
+       LaplaceDiffusionEntropyVariables3D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
        CompressibleNavierStokesDiffusion3D
 

src/equations/equations_parabolic.jl

@@ -14,6 +14,8 @@ include("laplace_diffusion_1d.jl")
 include("laplace_diffusion_2d.jl")
 include("laplace_diffusion_3d.jl")
 
+include("laplace_diffusion_entropy_variables.jl")
+
 # Compressible Navier-Stokes equations
 abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS, GradientVariables} <:
               AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} end

src/equations/laplace_diffusion_entropy_variables.jl

@@ -0,0 +1,121 @@
+@doc raw"""
+    LaplaceDiffusionEntropyVariables1D(equations)
+    LaplaceDiffusionEntropyVariables2D(equations)
+    LaplaceDiffusionEntropyVariables3D(equations)
+
+This represent an artificial viscosity term used to enforce entropy stability
+``\nabla \cdot (\epsilon(u)\frac{\partial u}{\partial w}\nabla w(u)))``, 
+where `w(u)` denotes the mapping between conservative and entropy variables. 
+"""
+struct LaplaceDiffusionEntropyVariables{NDIMS, E, N, T} <:
+       AbstractLaplaceDiffusion{NDIMS, N}
+    diffusivity::T
+    equations_hyperbolic::E
+end
+
+function LaplaceDiffusionEntropyVariables1D(diffusivity, equations_hyperbolic)
+    LaplaceDiffusionEntropyVariables{1, typeof(equations_hyperbolic),
+                                     nvariables(equations_hyperbolic),
+                                     typeof(diffusivity)}(diffusivity, equations_hyperbolic)
+end
+
+function LaplaceDiffusionEntropyVariables2D(diffusivity, equations_hyperbolic)
+    LaplaceDiffusionEntropyVariables{2, typeof(equations_hyperbolic),
+                                     nvariables(equations_hyperbolic),
+                                     typeof(diffusivity)}(diffusivity, equations_hyperbolic)
+end
+
+function LaplaceDiffusionEntropyVariables3D(diffusivity, equations_hyperbolic)
+    LaplaceDiffusionEntropyVariables{3, typeof(equations_hyperbolic),
+                                     nvariables(equations_hyperbolic),
+                                     typeof(diffusivity)}(diffusivity, equations_hyperbolic)
+end
+
+function varnames(variable_mapping, equations_parabolic::LaplaceDiffusionEntropyVariables)
+    varnames(variable_mapping, equations_parabolic.equations_hyperbolic)
+end
+
+function gradient_variable_transformation(::LaplaceDiffusionEntropyVariables)
+    cons2entropy
+end
+
+function cons2entropy(u, equations::LaplaceDiffusionEntropyVariables)
+    cons2entropy(u, equations.equations_hyperbolic)
+end
+function entropy2cons(w, equations::LaplaceDiffusionEntropyVariables)
+    entropy2cons(w, equations.equations_hyperbolic)
+end
+
+# generic fallback, assuming entropy2cons exists
+function jacobian_entropy2cons(w, equations)
+    return equations.diffusivity * ForwardDiff.jacobian(w -> entropy2cons(w, equations), w)
+end
+
+# Note that here, `u` should be the transformed entropy variables, and 
+# not the conservative variables.
+function flux(u, gradients, orientation::Integer,
+              equations::LaplaceDiffusionEntropyVariables{1})
+    dudx = gradients
+    diffusivity = jacobian_entropy2cons(u, equations)
+    # if orientation == 1
+    return SVector(diffusivity * dudx)
+end
+
+function flux(u, gradients, orientation::Integer,
+              equations::LaplaceDiffusionEntropyVariables{2})
+    dudx, dudy = gradients
+    diffusivity = jacobian_entropy2cons(u, equations)
+    if orientation == 1
+        return SVector(diffusivity * dudx)
+    else # if orientation == 2
+        return SVector(diffusivity * dudy)
+    end
+end
+
+function flux(u, gradients, orientation::Integer,
+              equations::LaplaceDiffusionEntropyVariables{3})
+    dudx, dudy, dudz = gradients
+    diffusivity = jacobian_entropy2cons(u, equations)
+    if orientation == 1
+        return SVector(diffusivity * dudx)
+    elseif orientation == 2
+        return SVector(diffusivity * dudy)
+    else # if orientation == 3
+        return SVector(diffusivity * dudz)
+    end
+end
+
+# Dirichlet and Neumann boundary conditions for use with parabolic solvers in weak form.
+# Note that these are general, so they apply to LaplaceDiffusionEntropyVariables in any 
+# spatial dimension.
+@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner,
+                                                                  normal::AbstractVector,
+                                                                  x, t,
+                                                                  operator_type::Gradient,
+                                                                  equations_parabolic::LaplaceDiffusionEntropyVariables)
+    return boundary_condition.boundary_value_function(x, t, equations_parabolic)
+end
+
+@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner,
+                                                                  normal::AbstractVector,
+                                                                  x, t,
+                                                                  operator_type::Divergence,
+                                                                  equations_parabolic::LaplaceDiffusionEntropyVariables)
+    return flux_inner
+end
+
+@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner,
+                                                                normal::AbstractVector,
+                                                                x, t,
+                                                                operator_type::Divergence,
+                                                                equations_parabolic::LaplaceDiffusionEntropyVariables)
+    return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic)
+end
+
+@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner,
+                                                                normal::AbstractVector,
+                                                                x, t,
+                                                                operator_type::Gradient,
+                                                                equations_parabolic::LaplaceDiffusionEntropyVariables)
+    return flux_inner
+end

test/test_dgmulti_2d.jl

@@ -940,6 +940,32 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_laplace_diffusion.jl (Tri, Polynomial)" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_laplace_diffusion.jl"),
+                        cells_per_dimension=8, element_type=Tri(),
+                        approximation_type=Polynomial(), tspan=(0.0, 0.2),
+                        l2=[
+                            0.09573222780014387,
+                            0.07305375992978504,
+                            0.07305375992978479,
+                            0.35291692218403836
+                        ],
+                        linf=[
+                            0.2632457353233918,
+                            0.22249956389208922,
+                            0.2224995638920848,
+                            0.9524896560641394
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 # Clean up afterwards: delete Trixi.jl output directory
 @test_nowarn isdir(outdir) && rm(outdir, recursive = true)

test/test_tree_1d_euler.jl

@@ -496,6 +496,26 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_euler_laplace_diffusion.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_laplace_diffusion.jl"),
+                        l2=[0.10954500481114468,
+                            0.1417583694046777,
+                            0.4087206508328759],
+                        linf=[
+                            0.17183237920520245,
+                            0.2203023610743297,
+                            0.6347464031934038
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 
 end # module

test/test_tree_3d_euler.jl

@@ -440,6 +440,32 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_laplace_diffusion.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_laplace_diffusion.jl"),
+                        l2=[
+                            0.013299230512542162,
+                            0.0073025819009651,
+                            0.007302581900965106,
+                            0.007300042097573285,
+                            0.04888085245959731
+                        ],
+                        linf=[
+                            0.31714843611640464,
+                            0.23586839231625517,
+                            0.23586839231625506,
+                            0.23698123351744782,
+                            1.1174271158464726
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
+
 @trixi_testset "elixir_euler_blob_amr.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_amr.jl"),
                         l2=[