Merge branch 'main' into sc/mhd-treemesh-3d

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/[email protected].23
+        uses: crate-ci/[email protected].26

NEWS.md

@@ -4,6 +4,12 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+## Changes in the v0.6 lifecycle
+
+#### Added
+- AMR for hyperbolic-parabolic equations on 3D `P4estMesh`
+- `flux_hllc` on non-cartesian meshes for `CompressibleEulerEquations{2,3}D`
+
 ## Changes when updating to v0.6 from v0.5.x
 
 #### Added
@@ -38,7 +44,7 @@ for human readability.
 - 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
-- Implementation of the quasi-1D shallow water equations
+- Implementation of the quasi-1D shallow water and compressible Euler equations
 - Subcell (positivity and local min/max) limiting support for conservative variables
   in 2D for `TreeMesh`
 - AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh`

Project.toml

@@ -1,11 +1,13 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.6.4-pre"
+version = "0.6.6-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
@@ -52,6 +54,8 @@ TrixiMakieExt = "Makie"
 [compat]
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
+DataStructures = "0.18.15"
+DiffEqBase = "6 - 6.143"
 DiffEqCallbacks = "2.25"
 EllipsisNotation = "1.0"
 FillArrays = "0.13.2, 1"
@@ -84,8 +88,8 @@ StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
 SummationByPartsOperators = "0.5.41"
-T8code = "0.4.3"
-TimerOutputs = "0.5"
+T8code = "0.4.3, 0.5"
+TimerOutputs = "0.5.7"
 Triangulate = "2.0"
 TriplotBase = "0.1"
 TriplotRecipes = "0.1"

examples/p4est_3d_dgsem/elixir_navierstokes_blast_wave_amr.jl

@@ -0,0 +1,113 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 6.25e-4 # equivalent to Re = 1600
+
+equations = CompressibleEulerEquations3D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+                                                          Prandtl = prandtl_number())
+
+function initial_condition_3d_blast_wave(x, t, equations::CompressibleEulerEquations3D)
+    rho_c = 1.0
+    p_c = 1.0
+    u_c = 0.0
+
+    rho_o = 0.125
+    p_o = 0.1
+    u_o = 0.0
+
+    rc = 0.5
+    r = sqrt(x[1]^2 + x[2]^2 + x[3]^2)
+    if r < rc
+        rho = rho_c
+        v1 = u_c
+        v2 = u_c
+        v3 = u_c
+        p = p_c
+    else
+        rho = rho_o
+        v1 = u_o
+        v2 = u_o
+        v3 = u_o
+        p = p_o
+    end
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+initial_condition = initial_condition_3d_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 1.0,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+coordinates_min = (-1.0, -1.0, -1.0) .* pi
+coordinates_max = (1.0, 1.0, 1.0) .* pi
+
+trees_per_dimension = (4, 4, 4)
+
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (true, true, true), initial_refinement_level = 1)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.8)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.05,
+                                      max_level = 3, max_threshold = 0.1)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 10,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        save_solution)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            ode_default_options()..., callback = callbacks)
+summary_callback() # print the timer summary

examples/p4est_3d_dgsem/elixir_navierstokes_taylor_green_vortex_amr.jl

@@ -0,0 +1,106 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Navier-Stokes equations
+
+# TODO: parabolic; unify names of these accessor functions
+prandtl_number() = 0.72
+mu() = 6.25e-4 # equivalent to Re = 1600
+
+equations = CompressibleEulerEquations3D(1.4)
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+                                                          Prandtl = prandtl_number())
+
+"""
+    initial_condition_taylor_green_vortex(x, t, equations::CompressibleEulerEquations3D)
+
+The classical Taylor-Green vortex.
+"""
+function initial_condition_taylor_green_vortex(x, t,
+                                               equations::CompressibleEulerEquations3D)
+    A = 1.0 # magnitude of speed
+    Ms = 0.1 # maximum Mach number
+
+    rho = 1.0
+    v1 = A * sin(x[1]) * cos(x[2]) * cos(x[3])
+    v2 = -A * cos(x[1]) * sin(x[2]) * cos(x[3])
+    v3 = 0.0
+    p = (A / Ms)^2 * rho / equations.gamma # scaling to get Ms
+    p = p +
+        1.0 / 16.0 * A^2 * rho *
+        (cos(2 * x[1]) * cos(2 * x[3]) + 2 * cos(2 * x[2]) + 2 * cos(2 * x[1]) +
+         cos(2 * x[2]) * cos(2 * x[3]))
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+initial_condition = initial_condition_taylor_green_vortex
+
+@inline function vel_mag(u, equations::CompressibleEulerEquations3D)
+    rho, rho_v1, rho_v2, rho_v3, _ = u
+    return sqrt(rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
+end
+
+volume_flux = flux_ranocha
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0, -1.0, -1.0) .* pi
+coordinates_max = (1.0, 1.0, 1.0) .* pi
+
+trees_per_dimension = (2, 2, 2)
+
+mesh = P4estMesh(trees_per_dimension, polydeg = 3,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = (true, true, true), initial_refinement_level = 0)
+
+semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
+                                             initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.5)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 50
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     save_analysis = true,
+                                     extra_analysis_integrals = (energy_kinetic,
+                                                                 energy_internal,
+                                                                 enstrophy))
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_indicator = IndicatorLöhner(semi, variable = vel_mag)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 0,
+                                      med_level = 1, med_threshold = 0.1,
+                                      max_level = 3, max_threshold = 0.2)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = false,
+                           adapt_initial_condition_only_refine = false)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        amr_callback,
+                        save_solution)
+
+###############################################################################
+# run the simulation
+
+time_int_tol = 1e-8
+sol = solve(ode, RDPK3SpFSAL49(); abstol = time_int_tol, reltol = time_int_tol,
+            ode_default_options()..., callback = callbacks)
+summary_callback() # print the timer summary

examples/tree_1d_dgsem/elixir_euler_quasi_1d_discontinuous.jl

@@ -0,0 +1,85 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# Semidiscretization of the quasi 1d compressible Euler equations
+# See Chan et al.  https://doi.org/10.48550/arXiv.2307.12089 for details
+
+equations = CompressibleEulerEquationsQuasi1D(1.4)
+
+"""
+    initial_condition_discontinuity(x, t, equations::CompressibleEulerEquations1D)
+
+A discontinuous initial condition taken from
+- Jesse Chan, Khemraj Shukla, Xinhui Wu, Ruofeng Liu, Prani Nalluri (2023)
+    High order entropy stable schemes for the quasi-one-dimensional
+    shallow water and compressible Euler equations
+    [DOI: 10.48550/arXiv.2307.12089](https://doi.org/10.48550/arXiv.2307.12089)   
+"""
+function initial_condition_discontinuity(x, t,
+                                         equations::CompressibleEulerEquationsQuasi1D)
+    rho = (x[1] < 0) ? 3.4718 : 2.0
+    v1 = (x[1] < 0) ? -2.5923 : -3.0
+    p = (x[1] < 0) ? 5.7118 : 2.639
+    a = (x[1] < 0) ? 1.0 : 1.5
+
+    return prim2cons(SVector(rho, v1, p, a), equations)
+end
+
+initial_condition = initial_condition_discontinuity
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_chan_etal)
+volume_flux = (flux_chan_etal, flux_nonconservative_chan_etal)
+
+basis = LobattoLegendreBasis(3)
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-1.0,)
+coordinates_max = (1.0,)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000)
+
+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,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+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