Merge branch 'main' into velocity

trixi-framework · Sep 30, 2024 · 82976ac · 82976ac
2 parents d29df2b + 1672c2a
commit 82976ac
Show file tree

Hide file tree

Showing 15 changed files with 185 additions and 40 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -11,6 +11,8 @@ for human readability.
 - 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/README.md b/README.md
@@ -37,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/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/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/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/src/callbacks_stage/subcell_limiter_idp_correction.jl b/src/callbacks_stage/subcell_limiter_idp_correction.jl
@@ -23,9 +23,6 @@ called with [`VolumeIntegralSubcellLimiting`](@ref).
 - Pazner (2020)
   Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting
   [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876)
-
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
 """
 struct SubcellLimiterIDPCorrection end
 

diff --git a/src/semidiscretization/semidiscretization_euler_gravity.jl b/src/semidiscretization/semidiscretization_euler_gravity.jl
@@ -134,7 +134,26 @@ function SemidiscretizationEulerGravity(semi_euler::SemiEuler,
                                                                       parameters, cache)
 end
 
-# TODO: AD, add appropriate method for remake
+function remake(semi::SemidiscretizationEulerGravity;
+                uEltype = real(semi.semi_gravity.solver),
+                semi_euler = semi.semi_euler,
+                semi_gravity = semi.semi_gravity,
+                parameters = semi.parameters)
+    semi_euler = remake(semi_euler, uEltype = uEltype)
+    semi_gravity = remake(semi_gravity, uEltype = uEltype)
+
+    # Recreate cache, i.e., registers for u with e.g. AD datatype
+    u_ode = compute_coefficients(zero(real(semi_gravity)), semi_gravity)
+    du_ode = similar(u_ode)
+    u_tmp1_ode = similar(u_ode)
+    u_tmp2_ode = similar(u_ode)
+    cache = (; u_ode, du_ode, u_tmp1_ode, u_tmp2_ode)
+
+    SemidiscretizationEulerGravity{typeof(semi_euler), typeof(semi_gravity),
+                                   typeof(parameters), typeof(cache)}(semi_euler,
+                                                                      semi_gravity,
+                                                                      parameters, cache)
+end
 
 function Base.show(io::IO, semi::SemidiscretizationEulerGravity)
     @nospecialize semi # reduce precompilation time

diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
@@ -196,9 +196,6 @@ with a low-order FV method. Used with limiter [`SubcellLimiterIDP`](@ref).
     mainly because the implementation assumes that low- and high-order schemes have the same
     surface terms, which is not guaranteed for non-conforming meshes. The low-order scheme
     with a high-order mortar is not invariant domain preserving.
-
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
 """
 struct VolumeIntegralSubcellLimiting{VolumeFluxDG, VolumeFluxFV, Limiter} <:
        AbstractVolumeIntegral

diff --git a/src/solvers/dgsem_tree/subcell_limiters.jl b/src/solvers/dgsem_tree/subcell_limiters.jl
@@ -60,9 +60,6 @@ where `d = #dimensions`). See equation (20) of Pazner (2020) and equation (30) o
 - Pazner (2020)
   Sparse invariant domain preserving discontinuous Galerkin methods with subcell convex limiting
   [DOI: 10.1016/j.cma.2021.113876](https://doi.org/10.1016/j.cma.2021.113876)
-
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
 """
 struct SubcellLimiterIDP{RealT <: Real, LimitingVariablesNonlinear,
                          LimitingOnesidedVariablesNonlinear, Cache} <:

diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl
@@ -19,9 +19,6 @@ The third-order SSP Runge-Kutta method of Shu and Osher.
 - Shu, Osher (1988)
   "Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes" (Eq. 2.18)
   [DOI: 10.1016/0021-9991(88)90177-5](https://doi.org/10.1016/0021-9991(88)90177-5)
-
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
 """
 struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP
     numerator_a::SVector{3, Float64}
@@ -133,9 +130,6 @@ end
 
 The following structures and methods provide the infrastructure for SSP Runge-Kutta methods
 of type `SimpleAlgorithmSSP`.
-
-!!! warning "Experimental implementation"
-    This is an experimental feature and may change in future releases.
 """
 function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP;
                dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)

diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl
@@ -68,15 +68,14 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages, eig_vals, tspan,
     a_matrix[:, 1] -= A
     a_matrix[:, 2] = A
 
-    return a_matrix, c
+    return a_matrix, c, dt_opt
 end
 
 # Compute the Butcher tableau for a paired explicit Runge-Kutta method order 2
 # using provided monomial coefficients file
 function compute_PairedExplicitRK2_butcher_tableau(num_stages,
                                                    base_path_monomial_coeffs::AbstractString,
                                                    bS, cS)
-
     # c Vector form Butcher Tableau (defines timestep per stage)
     c = zeros(num_stages)
     for k in 2:num_stages
@@ -107,7 +106,7 @@ function compute_PairedExplicitRK2_butcher_tableau(num_stages,
 end
 
 @doc raw"""
-    PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString,
+    PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString, dt_opt,
                       bS = 1.0, cS = 0.5)
     PairedExplicitRK2(num_stages, tspan, semi::AbstractSemidiscretization;
                       verbose = false, bS = 1.0, cS = 0.5)
@@ -118,6 +117,7 @@ end
     - `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing 
       monomial coefficients of the stability polynomial of PERK method.
       The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(num_stages).txt")` and separated by line breaks.
+    - `dt_opt` (`Float64`): Optimal time step size for the simulation setup.
     - `tspan`: Time span of the simulation.
     - `semi` (`AbstractSemidiscretization`): Semidiscretization setup.
     -  `eig_vals` (`Vector{ComplexF64}`): Eigenvalues of the Jacobian of the right-hand side (rhs) of the ODEProblem after the
@@ -144,16 +144,19 @@ mutable struct PairedExplicitRK2 <: AbstractPairedExplicitRKSingle
     b1::Float64
     bS::Float64
     cS::Float64
+    dt_opt::Float64
 end # struct PairedExplicitRK2
 
 # Constructor that reads the coefficients from a file
 function PairedExplicitRK2(num_stages, base_path_monomial_coeffs::AbstractString,
+                           dt_opt,
                            bS = 1.0, cS = 0.5)
+    # If the user has the monomial coefficients, they also must have the optimal time step
     a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages,
                                                             base_path_monomial_coeffs,
                                                             bS, cS)
 
-    return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS)
+    return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt)
 end
 
 # Constructor that calculates the coefficients with polynomial optimizer from a
@@ -171,12 +174,12 @@ end
 function PairedExplicitRK2(num_stages, tspan, eig_vals::Vector{ComplexF64};
                            verbose = false,
                            bS = 1.0, cS = 0.5)
-    a_matrix, c = compute_PairedExplicitRK2_butcher_tableau(num_stages,
-                                                            eig_vals, tspan,
-                                                            bS, cS;
-                                                            verbose)
+    a_matrix, c, dt_opt = compute_PairedExplicitRK2_butcher_tableau(num_stages,
+                                                                    eig_vals, tspan,
+                                                                    bS, cS;
+                                                                    verbose)
 
-    return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS)
+    return PairedExplicitRK2(num_stages, a_matrix, c, 1 - bS, bS, cS, dt_opt)
 end
 
 # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1
@@ -232,6 +235,26 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F,
     k_higher::uType
 end
 
+"""
+    calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode)
+
+This function computes the CFL number once using the initial condition of the problem and the optimal timestep (`dt_opt`) from the ODE algorithm.
+"""
+function calculate_cfl(ode_algorithm::AbstractPairedExplicitRKSingle, ode)
+    t0 = first(ode.tspan)
+    u_ode = ode.u0
+    semi = ode.p
+    dt_opt = ode_algorithm.dt_opt
+
+    mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
+    u = wrap_array(u_ode, mesh, equations, solver, cache)
+
+    cfl_number = dt_opt / max_dt(u, t0, mesh,
+                        have_constant_speed(equations), equations,
+                        solver, cache)
+    return cfl_number
+end
+
 """
     add_tstop!(integrator::PairedExplicitRK2Integrator, t)
 Add a time stop during the time integration process.

diff --git a/test/test_paper_self_gravitating_gas_dynamics.jl b/test/test_paper_self_gravitating_gas_dynamics.jl
@@ -215,16 +215,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_eulergravity_jeans_instability.jl"),
                         l2=[
-                            10733.63378538114,
-                            13356.780607423452,
+                            10733.28239179182,
+                            13356.0533511341,
                             1.6722844879795038e-6,
-                            26834.076821148774
+                            26833.19833691448
                         ],
                         linf=[
-                            15194.296424901113,
-                            18881.481685044182,
+                            15193.794080890715,
+                            18880.45819785685,
                             6.809726988008751e-6,
-                            37972.99700513482
+                            37971.74113135785
                         ],
                         tspan=(0.0, 0.1),
                         atol=4.0e-6)
@@ -242,16 +242,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_eulergravity_jeans_instability.jl"),
                         l2=[
-                            10734.598193238024,
-                            13358.217234481384,
+                            10734.653995035567,
+                            13357.709699808021,
                             1.911011743371934e-6,
-                            26836.487841241516
+                            26836.62734552835
                         ],
                         linf=[
-                            15195.661004798487,
-                            18883.512035906537,
+                            15195.73481107317,
+                            18882.799120551972,
                             7.867948710816926e-6,
-                            37976.408478975296
+                            37976.592992473394
                         ],
                         tspan=(0.0, 0.1),
                         atol=4.0e-6, # the background field is reatively large, so this corresponds to our usual atol

diff --git a/test/test_special_elixirs.jl b/test/test_special_elixirs.jl
@@ -254,6 +254,17 @@ end
                       tspan = (0.0, 0.0), initial_refinement_level = 0)
         @test_nowarn jacobian_ad_forward(semi)
     end
+
+    @timed_testset "EulerGravity" begin
+        trixi_include(@__MODULE__,
+                      joinpath(EXAMPLES_DIR,
+                               "paper_self_gravitating_gas_dynamics",
+                               "elixir_eulergravity_convergence.jl"),
+                      tspan = (0.0, 0.0), initial_refinement_level = 1)
+        J = jacobian_ad_forward(semi)
+        λ = eigvals(J)
+        @test maximum(real, λ) < 1.5
+    end
 end
 
 @timed_testset "Test linear structure (3D)" begin

diff --git a/test/test_tree_1d_advection.jl b/test/test_tree_1d_advection.jl
@@ -123,6 +123,25 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000
     end
 end
+
+# Testing the second-order paired explicit Runge-Kutta (PERK) method with the optimal CFL number
+@trixi_testset "elixir_advection_perk2_optimal_cfl.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_perk2_optimal_cfl.jl"),
+                        l2=[0.0009700887119146429],
+                        linf=[0.00137209242077041])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        # Larger values for allowed allocations due to usage of custom 
+        # integrator which are not *recorded* for the methods from 
+        # OrdinaryDiffEq.jl
+        # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 8000
+    end
+end
 end
 
 end # module