Skip to content

Commit aee3313

Browse files
jlchanranocha
jlchan
and
ranocha
authored
Add Laplace diffusion in terms of entropy variables (#2406)
* adding LaplaceDiffusionEntropyVariables type * removing export * adding Laplace entropy vars * rename and add diffusivity to equation type * add a test * format test * fix 3D equations * add 1D and 3D tests * epsilon -> kappa * separate files * formatting * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * add NEWS.md note * Update NEWS.md --------- Co-authored-by: Hendrik Ranocha <[email protected]>
1 parent a66a18f commit aee3313

13 files changed

+352
-0
lines changed

‎NEWS.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,8 @@ for human readability.
1010

1111
#### Added
1212

13+
- Added `LaplaceDiffusionEntropyVariables1D`, `LaplaceDiffusionEntropyVariables2D`, and `LaplaceDiffusionEntropyVariables3D`. These add scalar diffusion to each
14+
equation of a system, but apply diffusion in terms of the entropy variables, which symmetrizes the viscous formulation and ensures semi-discrete entropy dissipation ([#2406]).
1315
- Added the three-dimensional multi-ion magneto-hydrodynamics (MHD) equations with a
1416
generalized Lagrange multipliers (GLM) divergence cleaning technique ([#2215]).
1517
- New time integrator `PairedExplicitRK4`, implementing the fourth-order

‎examples/dgmulti_2d/elixir_euler_laplace_diffusion.jl‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
surface_flux = FluxLaxFriedrichs()
5+
volume_flux = flux_ranocha
6+
dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
7+
surface_integral = SurfaceIntegralWeakForm(surface_flux),
8+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
9+
10+
equations = CompressibleEulerEquations2D(1.4)
11+
equations_parabolic = LaplaceDiffusionEntropyVariables2D(0.001, equations)
12+
13+
initial_condition = initial_condition_weak_blast_wave
14+
15+
cells_per_dimension = (16, 16)
16+
mesh = DGMultiMesh(dg, cells_per_dimension,
17+
coordinates_min = (-1.0, -1.0), coordinates_max = (1.0, 1.0),
18+
periodicity = true)
19+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
20+
initial_condition, dg)
21+
22+
tspan = (0.0, 1.1)
23+
ode = semidiscretize(semi, tspan)
24+
25+
summary_callback = SummaryCallback()
26+
alive_callback = AliveCallback(alive_interval = 10)
27+
analysis_interval = 100
28+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg))
29+
callbacks = CallbackSet(summary_callback,
30+
analysis_callback,
31+
alive_callback)
32+
33+
###############################################################################
34+
# run the simulation
35+
36+
alg = RDPK3SpFSAL35()
37+
sol = solve(ode, alg; abstol = 1.0e-6, reltol = 1.0e-6,
38+
ode_default_options()..., callback = callbacks);

‎examples/tree_1d_dgsem/elixir_euler_laplace_diffusion.jl‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
###############################################################################
5+
# semidiscretization of the compressible Euler equations with Laplacian
6+
# diffusion represented in terms of entropy variables
7+
equations = CompressibleEulerEquations1D(1.4)
8+
equations_parabolic = LaplaceDiffusionEntropyVariables1D(0.01, equations)
9+
10+
initial_condition = initial_condition_weak_blast_wave
11+
12+
volume_flux = flux_ranocha
13+
solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
14+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
15+
16+
coordinates_min = (-2.0,)
17+
coordinates_max = (2.0,)
18+
mesh = TreeMesh(coordinates_min, coordinates_max,
19+
initial_refinement_level = 5,
20+
n_cells_max = 10_000)
21+
22+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
23+
initial_condition, solver)
24+
25+
###############################################################################
26+
# ODE solvers, callbacks etc.
27+
28+
tspan = (0.0, 0.4)
29+
ode = semidiscretize(semi, tspan)
30+
31+
summary_callback = SummaryCallback()
32+
33+
analysis_interval = 100
34+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
35+
36+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
37+
38+
save_solution = SaveSolutionCallback(interval = 100,
39+
save_initial_solution = true,
40+
save_final_solution = true,
41+
solution_variables = cons2prim)
42+
43+
stepsize_callback = StepsizeCallback(cfl = 0.8)
44+
45+
callbacks = CallbackSet(summary_callback,
46+
analysis_callback, alive_callback,
47+
save_solution,
48+
stepsize_callback)
49+
50+
###############################################################################
51+
# run the simulation
52+
53+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
54+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
55+
ode_default_options()..., callback = callbacks);

‎examples/tree_3d_dgsem/elixir_euler_laplace_diffusion.jl‎

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
###############################################################################
5+
# semidiscretization of the compressible Euler equations with Laplacian
6+
# diffusion represented in terms of entropy variables
7+
equations = CompressibleEulerEquations3D(1.4)
8+
equations_parabolic = LaplaceDiffusionEntropyVariables3D(0.01, equations)
9+
10+
initial_condition = initial_condition_weak_blast_wave
11+
12+
volume_flux = flux_ranocha
13+
solver = DGSEM(polydeg = 3, surface_flux = flux_ranocha,
14+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
15+
16+
coordinates_min = (-2.0, -2.0, -2.0)
17+
coordinates_max = (2.0, 2.0, 2.0)
18+
mesh = TreeMesh(coordinates_min, coordinates_max,
19+
initial_refinement_level = 3,
20+
n_cells_max = 100_000)
21+
22+
semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
23+
initial_condition, solver)
24+
25+
###############################################################################
26+
# ODE solvers, callbacks etc.
27+
28+
tspan = (0.0, 0.1)
29+
ode = semidiscretize(semi, tspan)
30+
31+
summary_callback = SummaryCallback()
32+
33+
analysis_interval = 100
34+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
35+
36+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
37+
38+
save_solution = SaveSolutionCallback(interval = 100,
39+
save_initial_solution = true,
40+
save_final_solution = true,
41+
solution_variables = cons2prim)
42+
43+
stepsize_callback = StepsizeCallback(cfl = 1.3)
44+
45+
callbacks = CallbackSet(summary_callback,
46+
analysis_callback, alive_callback,
47+
save_solution,
48+
stepsize_callback)
49+
50+
###############################################################################
51+
# run the simulation
52+
53+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
54+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
55+
ode_default_options()..., callback = callbacks);

‎src/Trixi.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -181,6 +181,8 @@ export AcousticPerturbationEquations2D,
181181
PassiveTracerEquations
182182

183183
export LaplaceDiffusion1D, LaplaceDiffusion2D, LaplaceDiffusion3D,
184+
LaplaceDiffusionEntropyVariables1D, LaplaceDiffusionEntropyVariables2D,
185+
LaplaceDiffusionEntropyVariables3D,
184186
CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
185187
CompressibleNavierStokesDiffusion3D
186188

‎src/equations/equations_parabolic.jl‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,11 @@ include("laplace_diffusion_1d.jl")
1414
include("laplace_diffusion_2d.jl")
1515
include("laplace_diffusion_3d.jl")
1616

17+
include("laplace_diffusion_entropy_variables.jl")
18+
include("laplace_diffusion_entropy_variables_1d.jl")
19+
include("laplace_diffusion_entropy_variables_2d.jl")
20+
include("laplace_diffusion_entropy_variables_3d.jl")
21+
1722
# Compressible Navier-Stokes equations
1823
abstract type AbstractCompressibleNavierStokesDiffusion{NDIMS, NVARS, GradientVariables} <:
1924
AbstractEquationsParabolic{NDIMS, NVARS, GradientVariables} end

‎src/equations/laplace_diffusion_entropy_variables.jl‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
@doc raw"""
2+
LaplaceDiffusionEntropyVariables1D(equations)
3+
LaplaceDiffusionEntropyVariables2D(equations)
4+
LaplaceDiffusionEntropyVariables3D(equations)
5+
6+
This represent a symmetrized Laplacian diffusion
7+
``\nabla \cdot (\kappa\frac{\partial u}{\partial w}\nabla w(u)))``,
8+
where ``w(u)`` denotes the mapping between conservative and entropy variables.
9+
Compared with `LaplaceDiffusion` (see [`LaplaceDiffusion1D`](@ref),
10+
[`LaplaceDiffusion2D`](@ref), and [`LaplaceDiffusion3D`](@ref)), `LaplaceDiffusionEntropyVariables` is
11+
guaranteed to dissipate entropy.
12+
"""
13+
struct LaplaceDiffusionEntropyVariables{NDIMS, E, N, T} <:
14+
AbstractLaplaceDiffusion{NDIMS, N}
15+
diffusivity::T
16+
equations_hyperbolic::E
17+
end
18+
19+
function varnames(variable_mapping, equations_parabolic::LaplaceDiffusionEntropyVariables)
20+
varnames(variable_mapping, equations_parabolic.equations_hyperbolic)
21+
end
22+
23+
function gradient_variable_transformation(::LaplaceDiffusionEntropyVariables)
24+
cons2entropy
25+
end
26+
27+
function cons2entropy(u, equations::LaplaceDiffusionEntropyVariables)
28+
cons2entropy(u, equations.equations_hyperbolic)
29+
end
30+
31+
function entropy2cons(w, equations::LaplaceDiffusionEntropyVariables)
32+
entropy2cons(w, equations.equations_hyperbolic)
33+
end
34+
35+
# This is used to compute the diffusivity tensor for LaplaceDiffusionEntropyVariables.
36+
# This is the generic fallback using AD (assuming entropy2cons exists)
37+
function jacobian_entropy2cons(w, equations)
38+
return equations.diffusivity * ForwardDiff.jacobian(w -> entropy2cons(w, equations), w)
39+
end
40+
41+
# Dirichlet and Neumann boundary conditions for use with parabolic solvers in weak form.
42+
# Note that these are general, so they apply to LaplaceDiffusionEntropyVariables in any
43+
# spatial dimension.
44+
@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner,
45+
normal::AbstractVector,
46+
x, t,
47+
operator_type::Gradient,
48+
equations_parabolic::LaplaceDiffusionEntropyVariables)
49+
return boundary_condition.boundary_value_function(x, t, equations_parabolic)
50+
end
51+
52+
@inline function (boundary_condition::BoundaryConditionDirichlet)(flux_inner, u_inner,
53+
normal::AbstractVector,
54+
x, t,
55+
operator_type::Divergence,
56+
equations_parabolic::LaplaceDiffusionEntropyVariables)
57+
return flux_inner
58+
end
59+
60+
@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner,
61+
normal::AbstractVector,
62+
x, t,
63+
operator_type::Divergence,
64+
equations_parabolic::LaplaceDiffusionEntropyVariables)
65+
return boundary_condition.boundary_normal_flux_function(x, t, equations_parabolic)
66+
end
67+
68+
@inline function (boundary_condition::BoundaryConditionNeumann)(flux_inner, u_inner,
69+
normal::AbstractVector,
70+
x, t,
71+
operator_type::Gradient,
72+
equations_parabolic::LaplaceDiffusionEntropyVariables)
73+
return flux_inner
74+
end

‎src/equations/laplace_diffusion_entropy_variables_1d.jl‎

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
function LaplaceDiffusionEntropyVariables1D(diffusivity, equations_hyperbolic)
2+
LaplaceDiffusionEntropyVariables{1, typeof(equations_hyperbolic),
3+
nvariables(equations_hyperbolic),
4+
typeof(diffusivity)}(diffusivity, equations_hyperbolic)
5+
end
6+
7+
# Note that here, `u` should be the transformed entropy variables, and
8+
# not the conservative variables.
9+
function flux(u, gradients, orientation::Integer,
10+
equations::LaplaceDiffusionEntropyVariables{1})
11+
dudx = gradients
12+
diffusivity = jacobian_entropy2cons(u, equations)
13+
# if orientation == 1
14+
return SVector(diffusivity * dudx)
15+
end

‎src/equations/laplace_diffusion_entropy_variables_2d.jl‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
function LaplaceDiffusionEntropyVariables2D(diffusivity, equations_hyperbolic)
2+
LaplaceDiffusionEntropyVariables{2, typeof(equations_hyperbolic),
3+
nvariables(equations_hyperbolic),
4+
typeof(diffusivity)}(diffusivity, equations_hyperbolic)
5+
end
6+
7+
function flux(u, gradients, orientation::Integer,
8+
equations::LaplaceDiffusionEntropyVariables{2})
9+
dudx, dudy = gradients
10+
diffusivity = jacobian_entropy2cons(u, equations)
11+
if orientation == 1
12+
return SVector(diffusivity * dudx)
13+
else # if orientation == 2
14+
return SVector(diffusivity * dudy)
15+
end
16+
end

‎src/equations/laplace_diffusion_entropy_variables_3d.jl‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
function LaplaceDiffusionEntropyVariables3D(diffusivity, equations_hyperbolic)
2+
LaplaceDiffusionEntropyVariables{3, typeof(equations_hyperbolic),
3+
nvariables(equations_hyperbolic),
4+
typeof(diffusivity)}(diffusivity, equations_hyperbolic)
5+
end
6+
7+
function flux(u, gradients, orientation::Integer,
8+
equations::LaplaceDiffusionEntropyVariables{3})
9+
dudx, dudy, dudz = gradients
10+
diffusivity = jacobian_entropy2cons(u, equations)
11+
if orientation == 1
12+
return SVector(diffusivity * dudx)
13+
elseif orientation == 2
14+
return SVector(diffusivity * dudy)
15+
else # if orientation == 3
16+
return SVector(diffusivity * dudz)
17+
end
18+
end

0 commit comments

Comments
 (0)