Sc/polytropic

NEWS.md

@@ -9,6 +9,7 @@ for human readability.
 #### Added
 
 - Experimental support for 3D parabolic diffusion terms has been added.
+- Implementation of compressible Euler equations in 2D
 
 #### Changed
 

examples/structured_2d_dgsem/elixir_euler_polytropic.jl

@@ -0,0 +1,88 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.001   # Isotropic simulation. gamma = 1.0 would lead to a singularity.
+# gamma = 2.0   # Adiabatic monatomic gas in 2d.
+kappa = 1.0     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+# Linear pressure wave in the negative x-direction.
+function initial_condition_wave(x, t, equations::PolytropicEulerEquations2D)
+  rho = 1.0
+  v1 = 0.0
+  if x[1] > 0.0
+    rho = ((1.0 + 0.01*sin(x[1]*2*pi)) / kappa)^(1/gamma)
+    v1 = ((0.01*sin((x[1]-1/2)*2*pi)) / kappa)
+  end
+  v2 = 0.0
+
+  return prim2cons(SVector(rho, v1, v2), equations)
+end
+initial_condition = initial_condition_wave
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg=2, surface_flux=flux_hll,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-2.0, -1.0)
+coordinates_max = ( 2.0,  1.0)
+
+cells_per_dimension = (64, 32)
+
+boundary_conditions = (
+                       x_neg=boundary_condition_periodic,
+                       x_pos=boundary_condition_periodic,
+                       y_neg=boundary_condition_periodic,
+                       y_pos=boundary_condition_periodic,
+                      )
+
+mesh = StructuredMesh(cells_per_dimension,
+                      coordinates_min,
+                      coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions=boundary_conditions)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 200
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval=analysis_interval)
+
+save_solution = SaveSolutionCallback(interval=50,
+                                     save_initial_solution=true,
+                                     save_final_solution=true,
+                                     solution_variables=cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl=1.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(1.0e-4, 1.0e-4),
+                                                     variables=(Trixi.density, pressure))
+
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition=false),
+            dt=1.0, # 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/structured_2d_dgsem/elixir_euler_polytropic_coupled.jl

@@ -0,0 +1,155 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Couple polytropic Euler systems.
+#
+# In a rectangular domain we solve two different sets of equations for the
+# left and the right half. Those are the isothermal equations (left)
+# and polytropic equations (right).
+# The coupling hapens on two interfaces. One is located at the center of the
+# domain such that the right half is coupled through its left boundary
+# and the left half is coupled through its right boundary. The other coupling
+# makes sure that the domain is periodic. Here the right boundary of the right
+# domain is coupled to the left boundary of the left domain.
+# The vertical boundaries are simply periodic.
+# As test case we send a linear wave through the domain and observe a change
+# in the dispersion relation when the wave enters the isothermal domain.
+
+###############################################################################
+# define the initial conditions
+
+function initial_condition_wave_isothermal(x, t, equations::PolytropicEulerEquations2D)
+    gamma = 1.001
+    kappa = 1.0
+
+    rho = 1.0
+    v1 = 0.0
+    v2 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2), equations)
+end
+
+function initial_condition_wave_polytropic(x, t, equations::PolytropicEulerEquations2D)
+    gamma = 2.0
+    kappa = 1.0
+
+    rho = 1.0
+    v1 = 0.0
+    if x[1] > 0.0
+        rho = ((1.0 + 0.01 * sin(x[1] * 2 * pi)) / kappa)^(1 / gamma)
+        v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / kappa)
+    end
+    v2 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2), equations)
+end
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+function Trixi.wrap_array(u_ode::AbstractVector, mesh::Trixi.AbstractMesh, equations,
+                          dg::DGSEM, cache)
+    Trixi.wrap_array_native(u_ode, mesh, equations, dg, cache)
+end
+
+# Define the equations involved.
+gamma1 = 1.001
+kappa1 = 1.0
+equations1 = PolytropicEulerEquations2D(gamma_A, kappa_A)
+gamma2 = 2.0
+kappa2 = 1.0
+equations2 = PolytropicEulerEquations2D(gamma_B, kappa_B)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg = 2, surface_flux = flux_hll,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min1 = (-2.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max1 = (0.0, 1.0) # maximum coordinates (max(x), max(y))
+coordinates_min2 = (0.0, -1.0) # minimum coordinates (min(x), min(y))
+coordinates_max2 = (2.0, 1.0) # maximum coordinates (max(x), max(y))
+
+cells_per_dimension = (32, 32)
+
+mesh1 = StructuredMesh(cells_per_dimension,
+                       coordinates_min1,
+                       coordinates_max1)
+mesh2 = StructuredMesh(cells_per_dimension,
+                       coordinates_min2,
+                       coordinates_max2)
+
+# A semidiscretization collects data structures and functions for the spatial discretization.
+semi1 = SemidiscretizationHyperbolic(mesh1, equations1,
+                                     initial_condition_wave_thermal, solver,
+                                     boundary_conditions = (x_neg = BoundaryConditionCoupled(2,
+                                                                                             (:end,
+                                                                                              :i_forward),
+                                                                                             Float64),
+                                                            x_pos = BoundaryConditionCoupled(2,
+                                                                                             (:begin,
+                                                                                              :i_forward),
+                                                                                             Float64),
+                                                            y_neg = boundary_condition_periodic,
+                                                            y_pos = boundary_condition_periodic))
+
+semi2 = SemidiscretizationHyperbolic(mesh2, equations2,
+                                     initial_condition_wave_polytropic, solver,
+                                     boundary_conditions = (x_neg = BoundaryConditionCoupled(1,
+                                                                                             (:end,
+                                                                                              :i_forward),
+                                                                                             Float64),
+                                                            x_pos = BoundaryConditionCoupled(1,
+                                                                                             (:begin,
+                                                                                              :i_forward),
+                                                                                             Float64),
+                                                            y_neg = boundary_condition_periodic,
+                                                            y_pos = boundary_condition_periodic))
+
+# Create a semidiscretization that bundles semi1 and semi2
+semi = SemidiscretizationCoupled(semi1, semi2)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 30.0
+ode = semidiscretize(semi, (0.0, 3.0))
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+resid_tol = 5.0e-12
+steady_state_callback = SteadyStateCallback(abstol = resid_tol, reltol = 0.0)
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback1 = AnalysisCallback(semi1, interval = 100)
+analysis_callback2 = AnalysisCallback(semi2, interval = 100)
+analysis_callback = AnalysisCallbackCoupled(semi, analysis_callback1, analysis_callback2)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval = 5,
+                                     solution_variables = cons2prim)
+
+# The StepsizeCallback handles the re-calculcation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+# 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)
+
+###############################################################################
+# 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 = 1.0, # 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()
+

src/Trixi.jl

@@ -146,7 +146,8 @@ export AcousticPerturbationEquations2D,
        LatticeBoltzmannEquations2D, LatticeBoltzmannEquations3D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
-       LinearizedEulerEquations2D
+       LinearizedEulerEquations2D,
+       PolytropicEulerEquations2D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D,
        CompressibleNavierStokesDiffusion2D, CompressibleNavierStokesDiffusion3D

src/equations/equations.jl

@@ -364,6 +364,11 @@ abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCO
 include("compressible_euler_multicomponent_1d.jl")
 include("compressible_euler_multicomponent_2d.jl")
 
+# PolytropicEulerEquations
+abstract type AbstractPolytropicEulerEquations{NDIMS, NVARS} <:
+              AbstractEquations{NDIMS, NVARS} end
+include("polytropic_euler_2d.jl")
+
 # Retrieve number of components from equation instance for the multicomponent case
 @inline function ncomponents(::AbstractCompressibleEulerMulticomponentEquations{NDIMS,
                                                                                 NVARS,