Merge branch 'main' into subcell-limiting-entropies

.github/workflows/Documenter.yml

@@ -34,7 +34,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1.9'
           show-versioninfo: true

.github/workflows/Downgrade.yml

@@ -64,7 +64,7 @@ jobs:
           - threaded
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

.github/workflows/Invalidations.yml

@@ -16,7 +16,7 @@ jobs:
     if: github.base_ref == github.event.repository.default_branch
     runs-on: ubuntu-latest
     steps:
-    - uses: julia-actions/setup-julia@v1
+    - uses: julia-actions/setup-julia@v2
       with:
         version: '1'
     - uses: actions/checkout@v4

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v4
       - name: Check spelling
-        uses: crate-ci/typos@v1.18.2
+        uses: crate-ci/typos@v1.19.0

.github/workflows/benchmark.yml

@@ -21,7 +21,7 @@ jobs:
       - run: |
           git fetch --tags
           git branch --create-reflog main origin/main
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

.github/workflows/ci.yml

@@ -107,7 +107,7 @@ jobs:
             trixi_test: threaded
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

.github/workflows/downstream.yml

@@ -64,7 +64,7 @@ jobs:
           - TrixiShallowWater.jl
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

NEWS.md

@@ -9,6 +9,8 @@ for human readability.
 #### Added
 - Implementation of `TimeSeriesCallback` for curvilinear meshes on `UnstructuredMesh2D` and extension
   to 1D and 3D on `TreeMesh`.
+- Implementation of 1D Linearized Euler Equations.
+- New analysis callback for 2D `P4estMesh` to compute integrated quantities along a boundary surface, e.g., pressure lift and drag coefficients.
 - Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh`.
 
 

Project.toml

@@ -1,7 +1,7 @@
 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.7.5-pre"
+version = "0.7.7-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
@@ -75,7 +75,7 @@ MuladdMacro = "0.2.2"
 Octavian = "0.3.21"
 OffsetArrays = "1.12"
 P4est = "0.4.9"
-Polyester = "0.7.5"
+Polyester = "0.7.10"
 PrecompileTools = "1.1"
 Preferences = "1.3"
 Printf = "1"

docs/src/performance.md

@@ -282,14 +282,3 @@ requires. It can thus be seen as a proxy for "energy used" and, as an extension,
     timing result, you need to set the analysis interval such that the
     `AnalysisCallback` is invoked at least once during the course of the simulation and
     discard the first PID value.
-
-## Performance issues with multi-threaded reductions
-[False sharing](https://en.wikipedia.org/wiki/False_sharing) is a known performance issue
-for systems with distributed caches. It also occurred for the implementation of a thread
-parallel bounds checking routine for the subcell IDP limiting
-in [PR #1736](https://github.com/trixi-framework/Trixi.jl/pull/1736).
-After some [testing and discussion](https://github.com/trixi-framework/Trixi.jl/pull/1736#discussion_r1423881895),
-it turned out that initializing a vector of length `n * Threads.nthreads()` and only using every
-n-th entry instead of a vector of length `Threads.nthreads()` fixes the problem.
-Since there are no processors with caches over 128B, we use `n = 128B / size(uEltype)`.
-Now, the bounds checking routine of the IDP limiting scales as hoped.

examples/dgmulti_2d/elixir_navierstokes_lid_driven_cavity.jl

@@ -4,12 +4,11 @@ using Trixi
 ###############################################################################
 # semidiscretization of the ideal compressible Navier-Stokes equations
 
-# TODO: parabolic; unify names of these accessor functions
 prandtl_number() = 0.72
-mu() = 0.001
+mu = 0.001
 
 equations = CompressibleEulerEquations2D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion2D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux

examples/dgmulti_3d/elixir_navierstokes_taylor_green_vortex.jl

@@ -5,12 +5,11 @@ 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
+mu = 6.25e-4 # equivalent to Re = 1600
 
 equations = CompressibleEulerEquations3D(1.4)
-equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu(),
+equations_parabolic = CompressibleNavierStokesDiffusion3D(equations, mu = mu,
                                                           Prandtl = prandtl_number())
 
 """

examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl

@@ -0,0 +1,132 @@
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+p_inf() = 1.0
+rho_inf() = p_inf() / (1.0 * 287.87) # p_inf = 1.0,  T = 1, R = 287.87
+mach_inf() = 0.85
+aoa() = pi / 180.0 # 1 Degree angle of attack
+c_inf(equations) = sqrt(equations.gamma * p_inf() / rho_inf())
+u_inf(equations) = mach_inf() * c_inf(equations)
+
+@inline function initial_condition_mach085_flow(x, t,
+                                                equations::CompressibleEulerEquations2D)
+    v1 = u_inf(equations) * cos(aoa())
+    v2 = u_inf(equations) * sin(aoa())
+
+    prim = SVector(rho_inf(), v1, v2, p_inf())
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_mach085_flow
+
+volume_flux = flux_ranocha_turbo
+surface_flux = flux_lax_friedrichs
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+shock_indicator = IndicatorHennemannGassner(equations, basis,
+                                            alpha_max = 0.5,
+                                            alpha_min = 0.001,
+                                            alpha_smooth = true,
+                                            variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(shock_indicator;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = volume_integral)
+
+mesh_file = Trixi.download("https://gist.githubusercontent.com/Arpit-Babbar/339662b4b46164a016e35c81c66383bb/raw/8bf94f5b426ba907ace87405cfcc1dcc2ef7cbda/NACA0012.inp",
+                           joinpath(@__DIR__, "NACA0012.inp"))
+
+mesh = P4estMesh{2}(mesh_file)
+
+# The outer boundary is constant but subsonic, so we cannot compute the
+# boundary flux for the external information alone. Thus, we use the numerical flux to distinguish
+# between inflow and outflow characteristics
+@inline function boundary_condition_subsonic_constant(u_inner,
+                                                      normal_direction::AbstractVector, x,
+                                                      t,
+                                                      surface_flux_function,
+                                                      equations::CompressibleEulerEquations2D)
+    u_boundary = initial_condition_mach085_flow(x, t, equations)
+
+    return Trixi.flux_hll(u_inner, u_boundary, normal_direction, equations)
+end
+
+boundary_conditions = Dict(:Left => boundary_condition_subsonic_constant,
+                           :Right => boundary_condition_subsonic_constant,
+                           :Top => boundary_condition_subsonic_constant,
+                           :Bottom => boundary_condition_subsonic_constant,
+                           :AirfoilBottom => boundary_condition_slip_wall,
+                           :AirfoilTop => boundary_condition_slip_wall)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers
+
+# Run for a long time to reach a steady state
+tspan = (0.0, 20.0)
+ode = semidiscretize(semi, tspan)
+
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 2000
+
+l_inf = 1.0 # Length of airfoil
+
+force_boundary_names = [:AirfoilBottom, :AirfoilTop]
+drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                           DragCoefficientPressure(aoa(), rho_inf(),
+                                                                   u_inf(equations), l_inf))
+
+lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+                                           LiftCoefficientPressure(aoa(), rho_inf(),
+                                                                   u_inf(equations), l_inf))
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     output_directory = "out",
+                                     save_analysis = true,
+                                     analysis_integrals = (drag_coefficient,
+                                                           lift_coefficient))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+amr_indicator = IndicatorLöhner(semi, variable = Trixi.density)
+
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 1,
+                                      med_level = 3, med_threshold = 0.05,
+                                      max_level = 4, max_threshold = 0.1)
+
+amr_interval = 100
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = amr_interval,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback, amr_callback)
+
+###############################################################################
+# run the simulation
+sol = solve(ode, SSPRK54(thread = OrdinaryDiffEq.True()),
+            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/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl

@@ -0,0 +1,125 @@
+using Downloads: download
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+@inline function initial_condition_mach038_flow(x, t,
+                                                equations::CompressibleEulerEquations2D)
+    # set the freestream flow parameters
+    rho_freestream = 1.4
+    v1 = 0.38
+    v2 = 0.0
+    p_freestream = 1.0
+
+    prim = SVector(rho_freestream, v1, v2, p_freestream)
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_mach038_flow
+
+volume_flux = flux_ranocha_turbo # FluxRotated(flux_chandrashekar) can also be used
+surface_flux = flux_lax_friedrichs
+
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+function mapping2cylinder(xi, eta)
+    xi_, eta_ = 0.5 * (xi + 1), 0.5 * (eta + 1.0) # Map from [-1,1] to [0,1] for simplicity
+
+    R2 = 50.0 # Bigger circle
+    R1 = 0.5  # Smaller circle
+
+    # Ensure an isotropic mesh by using elements with smaller radial length near the inner circle
+
+    r = R1 * exp(xi_ * log(R2 / R1))
+    theta = 2.0 * pi * eta_
+
+    x = r * cos(theta)
+    y = r * sin(theta)
+    return (x, y)
+end
+
+cells_per_dimension = (64, 64)
+# xi = -1 maps to the inner circle and xi = +1 maps to the outer circle and we can specify boundary
+# conditions there. However, the image of eta = -1, +1 coincides at the line y = 0. There is no
+# physical boundary there so we specify `periodicity = true` there and the solver treats the
+# element across eta = -1, +1 as neighbours which is what we want
+mesh = P4estMesh(cells_per_dimension, mapping = mapping2cylinder, polydeg = polydeg,
+                 periodicity = (false, true))
+
+# The boundary condition on the outer cylinder is constant but subsonic, so we cannot compute the
+# boundary flux from the external information alone. Thus, we use the numerical flux to distinguish
+# between inflow and outflow characteristics
+@inline function boundary_condition_subsonic_constant(u_inner,
+                                                      normal_direction::AbstractVector, x,
+                                                      t,
+                                                      surface_flux_function,
+                                                      equations::CompressibleEulerEquations2D)
+    u_boundary = initial_condition_mach038_flow(x, t, equations)
+
+    return surface_flux_function(u_inner, u_boundary, normal_direction, equations)
+end
+
+boundary_conditions = Dict(:x_neg => boundary_condition_slip_wall,
+                           :x_pos => boundary_condition_subsonic_constant)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers
+
+# Run for a long time to reach a steady state
+tspan = (0.0, 100.0)
+ode = semidiscretize(semi, tspan)
+
+# Callbacks
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 2000
+
+aoa = 0.0
+rho_inf = 1.4
+u_inf = 0.38
+l_inf = 1.0 # Diameter of circle
+
+drag_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
+                                           DragCoefficientPressure(aoa, rho_inf, u_inf,
+                                                                   l_inf))
+
+lift_coefficient = AnalysisSurfaceIntegral(semi, :x_neg,
+                                           LiftCoefficientPressure(aoa, rho_inf, u_inf,
+                                                                   l_inf))
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     output_directory = "out",
+                                     save_analysis = true,
+                                     analysis_integrals = (drag_coefficient,
+                                                           lift_coefficient))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 500,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 2.0)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+sol = solve(ode,
+            CarpenterKennedy2N54(williamson_condition = false;
+                                 thread = OrdinaryDiffEq.True()),
+            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