Sc/polytropic (#1526)

* Added polytropic Euler equations.

* Added polytropic Euler equations.

* Added polytropic Euler equations.

* Minor clean up in polytropic equations.

* Added example for the compressible polytropic Euler equations in 2d on a structured mesh.

* Cleaned up example elicir for polytropic Euler in 2d.

* Added polytrropic Euler in 2d in tests.

* Remoed unused routines.

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Added compressible Euler equations in news.

* Auto-reformatted equations.

* Auto-reformatted docs for polytropic equations.

* Added example elixir for polytropic and isothermal coupling.

* Corrected typos.

* Update examples/structured_2d_dgsem/elixir_euler_polytropic.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update NEWS.md

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Renamed elixir_euler_polytropic.jl to elixir_euler_polytropic.jl.

* Renamed Euler polytropic elixir in test.

* Renamed coupled polytropic example elixir.

* Reformatted the boundary conditions for the couple polytropic example elixir.

* Removed redundant inv_gamma_minus_one variable.

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl

Yes. That was necessary in an old version for Trixi.jl.

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_coupled.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Removed flux Ranocha for polytropic equations.

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Added Andrew's flux implementation flux_winters_etal.

* Added couped polytropic tests.

* Added Winter's fluxes for coupled polytropic elixir.

* Added flux_winters_etal.

* Added convergence test for polytropic Euler.

* Renamed polytropic uler convergence Elixir.

* Added convergence test initial condition.

* Corrected variable names in convergence test for polytropic equations.

* Update src/auxiliary/math.jl

Co-authored-by: Hendrik Ranocha <[email protected]>

* Added isotropic Eulre wave test case.

* Removed isotropic option in purely polytropic test case.

* Implmented Andrew's efficient way of computing coefficients.

* Added isotropic Euler wave as test case.

* Reformatted polytrropic_2d.

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Update src/auxiliary/math.jl

Co-authored-by: Andrew Winters <[email protected]>

* Reformatted example elixir for polytropic/isothermal equations.

* Updated errors for elixir_eulerpolytropic_coupled.jl.

* Added polytropic convergence test and renamed isothermal file.

* Corrected gamma and kappa.

* Added source terms for EOC.

* Added source term for ECO.

* Minor corrections.

* Added source terms.

* Corrected primary variables.

* Removed redundant kappa and gamme redefinitions.

* Improved style by using gamma and kappa from equations.

* Improved convergence test for polytropic.

* Corrected convergence test residual for polytropic equations.

* Allow use of flux_lax_friedrichs

* Fix comment

* Added flux calculation for polytropic Euler.

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_wave.jl

Co-authored-by: Andrew Winters <[email protected]>

* Moved initial condition and source term before flux function.

* Update examples/structured_2d_dgsem/elixir_eulerisothermal_wave.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Change kappa and gamma in some polytropic Euler elixirs.

* Correcte  source for Euler polytropic convergence test.

* Changed parameters gamma and kappa in Euler polytropic convergence test to other than 1.0.

* Reformatted polytropic Euler.

* Spelling errors.

* Reformatted.

* Added elixir for testing entropy conservation for polytropic Euler equations.

* Added weak blast wave as initial condition for the polytropic Euler equations.

* Added entropy conservation test.

* Typo in tests.

* Auto reformat.

* Corrected polytropic Euler test reults.

* Renamed entropy test for Euler polytropic.

* Renamed elixir for shock test.

* Added coupled Euler polytropic to the tests.

* Removed redundant flux function.

* Changed example parameters for coupled polytropic Euler.

* Removed max-speed calculations, as they are not being used.

* Removed coupled polytropic elixir, since this PR is on the polytropic
equations only.

* Correcte entropy variables for polytropic Euler.
This fixes the entropy conservation issue.

* Update examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Update src/equations/polytropic_euler_2d.jl

Co-authored-by: Andrew Winters <[email protected]>

* Format of comment.

* Update src/Trixi.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Updated polytropic wave test results after removing stage limiter.

* Removed stqage limiter.

* Reformatted Trixi.jl.

* Corrected removed exports in Trixi.jl.

* Fix for broadcasting.

* Renamed elixir for isothermal wave.

* Update test/test_structured_2d.jl

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

* Added comment about isothermal system in elixir.

---------

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>
Co-authored-by: Hendrik Ranocha <[email protected]>
Co-authored-by: Andrew Winters <[email protected]>

NEWS.md

@@ -12,6 +12,7 @@ for human readability.
 - Non-uniform `TreeMesh` available for hyperbolic-parabolic equations.
 - Capability to set truly discontinuous initial conditions in 1D.
 - Wetting and drying feature and examples for 1D and 2D shallow water equations
+- Implementation of the polytropic Euler equations in 2D
 - Subcell positivity limiting support for conservative variables in 2D for `TreeMesh`
 
 #### Changed

examples/structured_2d_dgsem/elixir_eulerpolytropic_convergence.jl

@@ -0,0 +1,63 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg = 3, surface_flux = flux_hll,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+
+cells_per_dimension = (4, 4)
+
+mesh = StructuredMesh(cells_per_dimension,
+                      coordinates_min,
+                      coordinates_max)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_convergence_test)
+
+###############################################################################
+# 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,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+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.1)
+
+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
+            save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_eulerpolytropic_ec.jl

@@ -0,0 +1,80 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.4
+kappa = 0.5     # Scaling factor for the pressure.
+equations = PolytropicEulerEquations2D(gamma, kappa)
+
+initial_condition = initial_condition_weak_blast_wave
+
+###############################################################################
+# Get the DG approximation space
+
+volume_flux = flux_winters_etal
+solver = DGSEM(polydeg=4, surface_flux=flux_winters_etal,
+               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))
+
+###############################################################################
+# Get the curved quad mesh from a mapping function
+
+# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D
+function mapping(xi_, eta_)
+  # Transform input variables between -1 and 1 onto [0,3]
+  xi = 1.5 * xi_ + 1.5
+  eta = 1.5 * eta_ + 1.5
+
+  y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) *
+                   cos(0.5 * pi * (2 * eta - 3)/3))
+
+  x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) *
+                  cos(2 * pi * (2 * y - 3)/3))
+
+  return SVector(x, y)
+end
+
+cells_per_dimension = (16, 16)
+
+# Create curved mesh with 16 x 16 elements
+mesh = StructuredMesh(cells_per_dimension, mapping)
+
+###############################################################################
+# create the semi discretization object
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+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)
+
+stepsize_callback = StepsizeCallback(cfl=1.0)
+
+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
+            save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_eulerpolytropic_isothermal_wave.jl

@@ -0,0 +1,83 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 1.0     # With gamma = 1 the system is isothermal.
+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)) / equations.kappa)^(1 / equations.gamma)
+        v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa)
+    end
+    v2 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2), equations)
+end
+initial_condition = initial_condition_wave
+
+volume_flux = flux_winters_etal
+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_eulerpolytropic_wave.jl

@@ -0,0 +1,80 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the polytropic Euler equations
+
+gamma = 2.0   # Adiabatic monatomic gas in 2d.
+kappa = 0.5   # 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)) / equations.kappa)^(1 / equations.gamma)
+        v1 = ((0.01 * sin((x[1] - 1 / 2) * 2 * pi)) / equations.kappa)
+    end
+    v2 = 0.0
+
+    return prim2cons(SVector(rho, v1, v2), equations)
+end
+initial_condition = initial_condition_wave
+
+volume_flux = flux_winters_etal
+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)
+
+###############################################################################
+# 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
+            save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()

src/Trixi.jl

@@ -150,7 +150,8 @@ export AcousticPerturbationEquations2D,
        ShallowWaterEquations1D, ShallowWaterEquations2D,
        ShallowWaterTwoLayerEquations1D, ShallowWaterTwoLayerEquations2D,
        ShallowWaterEquationsQuasi1D,
-       LinearizedEulerEquations2D
+       LinearizedEulerEquations2D,
+       PolytropicEulerEquations2D
 
 export LaplaceDiffusion1D, LaplaceDiffusion2D,
        CompressibleNavierStokesDiffusion1D, CompressibleNavierStokesDiffusion2D,
@@ -165,7 +166,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_kennedy_gruber, flux_shima_etal, flux_ec,
        flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal,
        flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,
-       flux_chan_etal, flux_nonconservative_chan_etal,
+       flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
 # TODO: TrixiShallowWater: move anything with "chen_noelle" to new file
        hydrostatic_reconstruction_chen_noelle, flux_nonconservative_chen_noelle,