Merge branch 'main' into surface_pressure_friction

trixi-framework · Oct 7, 2024 · 4745ec8 · 4745ec8
2 parents 9d7231d + 480aea9
commit 4745ec8
Show file tree

Hide file tree

Showing 127 changed files with 1,460 additions and 1,062 deletions.
diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml
@@ -29,7 +29,7 @@ jobs:
         # TODO: Change the call below to
         #       format(".")
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
           julia  -e 'using JuliaFormatter; format(["benchmark", "examples", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |

diff --git a/.github/workflows/SpellCheck.yml b/.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.24.3
+        uses: crate-ci/typos@v1.25.0
diff --git a/NEWS.md b/NEWS.md
@@ -4,13 +4,25 @@ 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 when updating to v0.9 from v0.8.x
+
+#### Changed
+
+- We removed the first argument `semi` corresponding to a `Semidiscretization` from the 
+  `AnalysisSurfaceIntegral` constructor, as it is no longer needed (see [#1959]).
+  The `AnalysisSurfaceIntegral` now only takes the arguments `boundary_symbols` and `variable`. ([#2069])
+- In functions `rhs!`, `rhs_parabolic!`  we removed the unused argument `initial_condition`. ([#2037])
+  Users should not be affected by this.
+
 ## Changes in the v0.8 lifecycle
 
 #### Changed
 
 - The AMR routines for `P4estMesh` and `T8codeMesh` were changed to work on the product
   of the Jacobian and the conserved variables instead of the conserved variables only
   to make AMR fully conservative ([#2028]). This may change AMR results slightly.
+- Subcell (IDP) limiting is now officially supported and not marked as experimental
+  anymore (see `VolumeIntegralSubcellLimiting`).
 
 ## Changes when updating to v0.8 from v0.7.x
 

diff --git a/Project.toml b/Project.toml
@@ -1,9 +1,10 @@
 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.8.9-DEV"
+version = "0.8.11-DEV"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -60,6 +61,7 @@ TrixiConvexECOSExt = ["Convex", "ECOS"]
 TrixiMakieExt = "Makie"
 
 [compat]
+Accessors = "0.1.12"
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
 Convex = "0.16"

diff --git a/README.md b/README.md
@@ -19,28 +19,6 @@
   <img width="300px" src="https://trixi-framework.github.io/assets/logo.png">
 </p>
 
-***
-**Trixi.jl at JuliaCon 2024**<br/>
-At this year's JuliaCon in Eindhoven, Netherlands, we will be present with several contributions
-from the Trixi Framework ecosystem:
-
-* [**Julia for Particle-Based Multiphysics with TrixiParticles.jl**](https://pretalx.com/juliacon2024/talk/TPFF8L/),<br/>
-  [*Erik Faulhaber*](https://github.com/efaulhaber/), [*Niklas Neher*](https://github.com/lasnikas/),
-  10th July 2024, 11:00–11:30, Function (4.1)
-* [**Towards Aerodynamic Simulations in Julia with Trixi.jl**](https://pretalx.com/juliacon2024/talk/XH8KBG/),<br/>
-  [*Daniel Doehring*](https://github.com/danieldoehring/),
-  10th July 2024, 15:30–16:00, While Loop (4.2)
-* [**libtrixi: serving legacy codes in earth system modeling with fresh Julia CFD**](https://pretalx.com/juliacon2024/talk/SXC7LA/),<br/>
-  [*Benedict Geihe*](https://github.com/benegee/),
-  12th July 2024, 14:00–17:00, Function (4.1)
-
-The last talk is part of the 
-[Julia for High-Performance Computing](https://juliacon.org/2024/minisymposia/hpc/)
-minisymposium, which this year is hosted by our own [*Hendrik Ranocha*](https://github.com/ranocha/).
-
-We are looking forward to seeing you there ♥️
-***
-
 **Trixi.jl** is a numerical simulation framework for conservation
 laws written in [Julia](https://julialang.org). A key objective for the
 framework is to be useful to both scientists and students. Therefore, next to
@@ -59,6 +37,7 @@ installation and postprocessing procedures. Its features include:
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
   * Positivity-preserving limiting
+  * Subcell invariant domain-preserving (IDP) limiting
   * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl)
 * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/)
   * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods)

diff --git a/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl b/docs/literate/src/files/behind_the_scenes_simulation_setup_plots/src/rhs_structure_figure.jl
@@ -35,7 +35,7 @@ plot!(Shape([(-1.7,-11), (1.7,-11), (1.7,-17), (-1.7,-17)]), linecolor=:blue, fi
 annotate!(1.5, -14, ("Trixi.jl", 12, :blue, :center))
 
 plot!(Shape([(-1.4,-14.5), (1.4,-14.5), (1.4,-16.5), (-1.4,-16.5)]), linecolor="black", fillcolor="white", label=false,linewidth=2)
-annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, initial_condition, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
+annotate!(0, -15.5, ("Trixi.rhs!(du, u, t, mesh, equations, \nboundary_conditions, source_terms, dg, cache)", 12, :black, :center))
 
 
 

diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl
@@ -51,7 +51,11 @@
 # The Newton-bisection algorithm is an iterative method and requires some parameters.
 # It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and
 # relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are
-# sufficient in most cases and therefore used as default. Additionally, there is the parameter
+# sufficient in most cases and therefore used as default. If the implemented bounds checking
+# functionality indicates problems with the limiting (see [below](@ref subcell_bounds_check))
+# the Newton method with the chosen parameters might not manage to converge. If so, adapting
+# the mentioned parameters helps fix that.
+# Additionally, there is the parameter
 # `gamma_constant_newton`, which can be used to scale the antidiffusive flux for the computation
 # of the blending coefficients of nonlinear variables. The default value is `2 * ndims(equations)`,
 # as it was shown by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) [Section 4.2.2.]
@@ -244,7 +248,7 @@ plot(sol)
 # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116)
 
 
-# ## Bounds checking
+# ## [Bounds checking](@id subcell_bounds_check)
 # Subcell limiting is based on the fulfillment of target bounds - either global or local.
 # Although the implementation works and has been thoroughly tested, there are some cases where
 # these bounds are not met.

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -30,6 +30,7 @@ installation and postprocessing procedures. Its features include:
   * Kinetic energy-preserving and entropy-stable methods based on flux differencing
   * Entropy-stable shock capturing
   * Positivity-preserving limiting
+  * Subcell invariant domain-preserving (IDP) limiting
   * [Finite difference summation by parts (SBP) methods](https://github.com/ranocha/SummationByPartsOperators.jl)
 * Compatible with the [SciML ecosystem for ordinary differential equations](https://diffeq.sciml.ai/latest/)
   * [Explicit low-storage Runge-Kutta time integration](https://diffeq.sciml.ai/latest/solvers/ode_solve/#Low-Storage-Methods)

diff --git a/docs/src/styleguide.md b/docs/src/styleguide.md
@@ -51,7 +51,7 @@ PRs that verify that running JuliaFormatter.jl again will not change the source
 To format your contributions before created a PR (or, at least, before requesting a review
 of your PR), you need to install JuliaFormatter.jl first by running
 ```shell
-julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
 ```
 You can then recursively format the core Julia files in the Trixi.jl repo by executing
 ```shell

diff --git a/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl b/examples/p4est_2d_dgsem/elixir_euler_NACA0012airfoil_mach085.jl
@@ -85,11 +85,11 @@ analysis_interval = 2000
 l_inf = 1.0 # Length of airfoil
 
 force_boundary_names = (:AirfoilBottom, :AirfoilTop)
-drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations), l_inf))
 

diff --git a/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl b/examples/p4est_2d_dgsem/elixir_euler_subsonic_cylinder.jl
@@ -89,11 +89,11 @@ rho_inf = 1.4
 u_inf = 0.38
 l_inf = 1.0 # Diameter of circle
 
-drag_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
+drag_coefficient = AnalysisSurfaceIntegral((:x_neg,),
                                            DragCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, (:x_neg,),
+lift_coefficient = AnalysisSurfaceIntegral((:x_neg,),
                                            LiftCoefficientPressure(aoa, rho_inf, u_inf,
                                                                    l_inf))
 

diff --git a/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl b/examples/p4est_2d_dgsem/elixir_navierstokes_NACA0012airfoil_mach08.jl
@@ -120,23 +120,23 @@ summary_callback = SummaryCallback()
 analysis_interval = 2000
 
 force_boundary_names = (:AirfoilBottom, :AirfoilTop)
-drag_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            DragCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),
                                                                    l_inf()))
 
-lift_coefficient = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient = AnalysisSurfaceIntegral(force_boundary_names,
                                            LiftCoefficientPressure(aoa(), rho_inf(),
                                                                    u_inf(equations),
                                                                    l_inf()))
 
-drag_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+drag_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                        DragCoefficientShearStress(aoa(),
                                                                                   rho_inf(),
                                                                                   u_inf(equations),
                                                                                   l_inf()))
 
-lift_coefficient_shear_force = AnalysisSurfaceIntegral(semi, force_boundary_names,
+lift_coefficient_shear_force = AnalysisSurfaceIntegral(force_boundary_names,
                                                        LiftCoefficientShearStress(aoa(),
                                                                                   rho_inf(),
                                                                                   u_inf(equations),

diff --git a/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl b/examples/paper_self_gravitating_gas_dynamics/elixir_eulergravity_jeans_instability.jl
@@ -91,7 +91,7 @@ semi_gravity = SemidiscretizationHyperbolic(mesh, equations_gravity, initial_con
 # combining both semidiscretizations for Euler + self-gravity
 parameters = ParametersEulerGravity(background_density = 1.5e7, # aka rho0
                                     gravitational_constant = 6.674e-8, # aka G
-                                    cfl = 1.6,
+                                    cfl = 0.8, # value as used in the paper
                                     resid_tol = 1.0e-4,
                                     n_iterations_max = 1000,
                                     timestep_gravity = timestep_gravity_carpenter_kennedy_erk54_2N!)
@@ -105,7 +105,7 @@ ode = semidiscretize(semi, tspan);
 
 summary_callback = SummaryCallback()
 
-stepsize_callback = StepsizeCallback(cfl = 1.0)
+stepsize_callback = StepsizeCallback(cfl = 0.5) # value as used in the paper
 
 save_solution = SaveSolutionCallback(interval = 10,
                                      save_initial_solution = true,

diff --git a/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl
@@ -51,7 +51,9 @@ limiter_idp = SubcellLimiterIDP(equations, basis;
                                 local_twosided_variables_cons = ["rho"],
                                 local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
                                                                        min)],
-                                max_iterations_newton = 40, # Default value of 10 iterations is too low to fulfill bounds.
+                                # Default parameters are not sufficient to fulfill bounds properly.
+                                max_iterations_newton = 40,
+                                newton_tolerances = (1.0e-13, 1.0e-15),
                                 positivity_variables_cons = [],
                                 positivity_variables_nonlinear = [])
 # Variables for global limiting (`positivity_variables_cons` and
@@ -96,7 +98,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.7)
+stepsize_callback = StepsizeCallback(cfl = 0.6)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

diff --git a/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl b/examples/tree_1d_dgsem/elixir_advection_perk2_optimal_cfl.jl
@@ -0,0 +1,84 @@
+
+using Convex, ECOS
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the linear advection equation
+
+advection_velocity = 1.0
+equations = LinearScalarAdvectionEquation1D(advection_velocity)
+
+# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
+solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
+
+coordinates_min = -1.0 # minimum coordinate
+coordinates_max = 1.0 # maximum coordinate
+
+# Create a uniformly refined mesh with periodic boundaries
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 30_000) # set maximum capacity of tree data structure
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,
+                                    solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 20.0
+tspan = (0.0, 20.0)
+ode = semidiscretize(semi, tspan);
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(alive_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 0.1,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 4,
+                                      med_level = 5, med_threshold = 0.1,
+                                      max_level = 6, max_threshold = 0.6)
+
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+# Construct second order paired explicit Runge-Kutta method with 6 stages for given simulation setup.
+# Pass `tspan` to calculate maximum time step allowed for the bisection algorithm used 
+# in calculating the polynomial coefficients in the ODE algorithm.
+ode_algorithm = Trixi.PairedExplicitRK2(6, tspan, semi)
+
+# For Paired Explicit Runge-Kutta methods, we receive an optimized timestep for a given reference semidiscretization.
+# To allow for e.g. adaptivity, we reverse-engineer the corresponding CFL number to make it available during the simulation.
+cfl_number = Trixi.calculate_cfl(ode_algorithm, ode)
+stepsize_callback = StepsizeCallback(cfl = cfl_number)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback,
+                        alive_callback,
+                        save_solution,
+                        analysis_callback,
+                        amr_callback,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+sol = Trixi.solve(ode, ode_algorithm,
+                  dt = 1.0, # Manual time step value, will be overwritten by the stepsize_callback when it is specified.
+                  save_everystep = false, callback = callbacks);
+
+# Print the timer summary
+summary_callback()
diff --git a/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl b/examples/tree_1d_dgsem/elixir_burgers_rarefaction.jl
@@ -43,10 +43,7 @@ end
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-function inflow(x, t, equations::InviscidBurgersEquation1D)
-    return initial_condition_rarefaction(coordinate_min, t, equations)
-end
-boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_rarefaction)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,

diff --git a/examples/tree_1d_dgsem/elixir_burgers_shock.jl b/examples/tree_1d_dgsem/elixir_burgers_shock.jl
@@ -43,10 +43,7 @@ end
 ###############################################################################
 # Specify non-periodic boundary conditions
 
-function inflow(x, t, equations::InviscidBurgersEquation1D)
-    return initial_condition_shock(coordinate_min, t, equations)
-end
-boundary_condition_inflow = BoundaryConditionDirichlet(inflow)
+boundary_condition_inflow = BoundaryConditionDirichlet(initial_condition_shock)
 
 function boundary_condition_outflow(u_inner, orientation, normal_direction, x, t,
                                     surface_flux_function,

diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl
@@ -43,7 +43,10 @@ basis = LobattoLegendreBasis(3)
 limiter_idp = SubcellLimiterIDP(equations, basis;
                                 local_twosided_variables_cons = ["rho"],
                                 local_onesided_variables_nonlinear = [(Trixi.entropy_math,
-                                                                       max)])
+                                                                       max)],
+                                # Default parameters are not sufficient to fulfill bounds properly.
+                                max_iterations_newton = 70,
+                                newton_tolerances = (1.0e-13, 1.0e-14))
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)