diff --git a/examples/tree_1d_dgsem/elixir_euler_adaptive_volume_integral.jl b/examples/tree_1d_dgsem/elixir_euler_adaptive_volume_integral.jl new file mode 100644 index 00000000000..44893a05b0f --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_adaptive_volume_integral.jl @@ -0,0 +1,58 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerEquations1D(1.4) + +initial_condition = initial_condition_weak_blast_wave + +basis = LobattoLegendreBasis(5) +surface_flux = flux_lax_friedrichs +# Note: Plain weak form volume integral is still stable for this problem + +volume_integral_fluxdiff = VolumeIntegralFluxDifferencing(flux_ranocha) + +# Standard Flux-Differencing volume integral, roughly twice as expensive as the adaptive one +#solver = DGSEM(basis, surface_flux, volume_integral_fluxdiff) + +indicator = IndicatorEntropyViolation(basis; threshold = 1e-6) +volume_integral = VolumeIntegralAdaptive(indicator; + volume_integral_default = VolumeIntegralWeakForm(), + volume_integral_stabilized = volume_integral_fluxdiff) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-2.0,) +coordinates_max = (2.0,) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 8, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.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) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false); + dt = 1.0, + ode_default_options()..., callback = callbacks); diff --git a/examples/tree_2d_dgsem/elixir_euler_warm_bubble_adaptive_integral.jl b/examples/tree_2d_dgsem/elixir_euler_warm_bubble_adaptive_integral.jl new file mode 100644 index 00000000000..8c16a3cb97d --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_euler_warm_bubble_adaptive_integral.jl @@ -0,0 +1,148 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +# Warm bubble test case from +# - Wicker, L. J., and Skamarock, W. C. (1998) +# A time-splitting scheme for the elastic equations incorporating +# second-order Runge–Kutta time differencing +# [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2) +# See also +# - Bryan and Fritsch (2002) +# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models +# [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2) +# - Carpenter, Droegemeier, Woodward, Hane (1990) +# Application of the Piecewise Parabolic Method (PPM) to +# Meteorological Modeling +# [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2) +struct WarmBubbleSetup + # Physical constants + g::Float64 # gravity of earth + c_p::Float64 # heat capacity for constant pressure (dry air) + c_v::Float64 # heat capacity for constant volume (dry air) + gamma::Float64 # heat capacity ratio (dry air) + + function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v) + new(g, c_p, c_v, gamma) + end +end + +# Initial condition +function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquations2D) + RealT = eltype(x) + @unpack g, c_p, c_v = setup + + # center of perturbation + center_x = 10000 + center_z = 2000 + # radius of perturbation + radius = 2000 + # distance of current x to center of perturbation + r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2) + + # perturbation in potential temperature + potential_temperature_ref = 300 + potential_temperature_perturbation = zero(RealT) + if r <= radius + potential_temperature_perturbation = 2 * cospi(0.5f0 * r / radius)^2 + end + potential_temperature = potential_temperature_ref + potential_temperature_perturbation + + # Exner pressure, solves hydrostatic equation for x[2] + exner = 1 - g / (c_p * potential_temperature) * x[2] + + # pressure + p_0 = 100_000 # reference pressure + R = c_p - c_v # gas constant (dry air) + p = p_0 * exner^(c_p / R) + + # temperature + T = potential_temperature * exner + + # density + rho = p / (R * T) + + v1 = 20 + v2 = 0 + E = c_v * T + 0.5f0 * (v1^2 + v2^2) + return SVector(rho, rho * v1, rho * v2, rho * E) +end + +# Source terms +@inline function (setup::WarmBubbleSetup)(u, x, t, equations::CompressibleEulerEquations2D) + @unpack g = setup + rho, _, rho_v2, _ = u + return SVector(zero(eltype(u)), zero(eltype(u)), -g * rho, -g * rho_v2) +end + +############################################################################### +# semidiscretization of the compressible Euler equations +warm_bubble_setup = WarmBubbleSetup() + +equations = CompressibleEulerEquations2D(warm_bubble_setup.gamma) + +boundary_conditions = (x_neg = boundary_condition_periodic, + x_pos = boundary_condition_periodic, + y_neg = boundary_condition_slip_wall, + y_pos = boundary_condition_slip_wall) + +polydeg = 5 +basis = LobattoLegendreBasis(polydeg) + +# This is a good estimate for the speed of sound in this example. +# Other values between 300 and 400 should work as well. +surface_flux = FluxLMARS(340.0) + +volume_flux = flux_kennedy_gruber + +indicator = IndicatorEntropyViolation(basis; threshold = 1e-3) +volume_integral = VolumeIntegralAdaptive(indicator; + volume_integral_default = VolumeIntegralWeakForm(), + volume_integral_stabilized = VolumeIntegralFluxDifferencing(volume_flux)) + +# This would be the standard version +#volume_integral = VolumeIntegralFluxDifferencing(volume_flux) + +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, -5000.0) +coordinates_max = (20_000.0, 15_000.0) + +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10_000, + periodicity = (true, false)) + +semi = SemidiscretizationHyperbolic(mesh, equations, warm_bubble_setup, solver, + source_terms = warm_bubble_setup, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1000.0) # seconds + +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 10_000 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = analysis_interval, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, + alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false, thread = Trixi.True()); + dt = 1.0, # solve needs a value here, will be overwritten by the stepsize_callback + ode_default_options()..., callback = callbacks); diff --git a/src/Trixi.jl b/src/Trixi.jl index 6b596bd74fb..be41ebc3bbd 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -262,6 +262,7 @@ export DG, VolumeIntegralFluxDifferencing, VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2, VolumeIntegralShockCapturingHG, IndicatorHennemannGassner, + VolumeIntegralAdaptive, IndicatorEntropyViolation, VolumeIntegralUpwind, SurfaceIntegralWeakForm, SurfaceIntegralStrongForm, SurfaceIntegralUpwind, diff --git a/src/basic_types.jl b/src/basic_types.jl index 1aaf6d6d464..7c4b2cd9315 100644 --- a/src/basic_types.jl +++ b/src/basic_types.jl @@ -49,6 +49,10 @@ abstract type AdaptorAMR{RealT <: Real} end # which will be specialized for different SBP bases abstract type AdaptorL2{RealT <: Real} <: AdaptorAMR{RealT} end +# abstract supertype of indicators used for AMR and/or shock capturing via +# volume-integral selection +abstract type AbstractIndicator end + # TODO: Taal decide, which abstract types shall be defined here? struct BoundaryConditionPeriodic end diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl index d2649384a18..27927e8dbb3 100644 --- a/src/solvers/dg.jl +++ b/src/solvers/dg.jl @@ -177,8 +177,60 @@ function Base.show(io::IO, mime::MIME"text/plain", end function get_element_variables!(element_variables, u, mesh, equations, - volume_integral::VolumeIntegralShockCapturingHG, dg, - cache) + volume_integral::VolumeIntegralShockCapturingHG, + dg, cache) + # call the indicator to get up-to-date values for IO + volume_integral.indicator(u, mesh, equations, dg, cache) + get_element_variables!(element_variables, volume_integral.indicator, + volume_integral) +end + +""" + VolumeIntegralAdaptive(indicator; + volume_integral_default = VolumeIntegralWeakForm(), + volume_integral_stabilized = VolumeIntegralFluxDifferencing(flux_central)) + +Possible combination: [`VolumeIntegralWeakForm`](@ref) and [`VolumeIntegralFluxDifferencing`](@ref). +TODO: Extend to [`VolumeIntegralWeakForm`](@ref) and [`VolumeIntegralShockCapturingHG`](@ref). +""" +struct VolumeIntegralAdaptive{VolumeIntegralDefault, VolumeIntegralStabilized, + Indicator} <: AbstractVolumeIntegral + volume_integral_default::VolumeIntegralDefault # Cheap(er) default volume integral to be used in non-critical regions + volume_integral_stabilized::VolumeIntegralStabilized # More expensive volume integral with stabilizing effect + indicator::Indicator +end + +function VolumeIntegralAdaptive(indicator; + volume_integral_default = VolumeIntegralWeakForm(), + volume_integral_stabilized = VolumeIntegralFluxDifferencing(flux_central)) + VolumeIntegralAdaptive{typeof(volume_integral_default), + typeof(volume_integral_stabilized), + typeof(indicator)}(volume_integral_default, + volume_integral_stabilized, + indicator) +end + +function Base.show(io::IO, mime::MIME"text/plain", + integral::VolumeIntegralAdaptive) + @nospecialize integral # reduce precompilation time + + if get(io, :compact, false) + show(io, integral) + else + summary_header(io, "VolumeIntegralAdaptive") + summary_line(io, "volume integral default", + integral.volume_integral_default) + summary_line(io, "volume integral stabilized", + integral.volume_integral_stabilized) + summary_line(io, "indicator", integral.indicator |> typeof |> nameof) + show(increment_indent(io), mime, integral.indicator) + summary_footer(io) + end +end + +function get_element_variables!(element_variables, u, mesh, equations, + volume_integral::VolumeIntegralAdaptive, + dg, cache) # call the indicator to get up-to-date values for IO volume_integral.indicator(u, mesh, equations, dg, cache) get_element_variables!(element_variables, volume_integral.indicator, diff --git a/src/solvers/dgsem/calc_volume_integral.jl b/src/solvers/dgsem/calc_volume_integral.jl index 5739a28b202..89200bb8bd6 100644 --- a/src/solvers/dgsem/calc_volume_integral.jl +++ b/src/solvers/dgsem/calc_volume_integral.jl @@ -11,6 +11,13 @@ function create_cache(mesh, equations, return NamedTuple() end +function create_cache(mesh, equations, + volume_integral::VolumeIntegralAdaptive, dg::DG, uEltype) + return create_cache(mesh, equations, + volume_integral.volume_integral_stabilized, + dg, uEltype) +end + # The following `calc_volume_integral!` functions are # dimension and meshtype agnostic, i.e., valid for all 1D, 2D, and 3D meshes. @@ -79,6 +86,41 @@ function calc_volume_integral!(du, u, mesh, return nothing end +function calc_volume_integral!(du, u, mesh, + nonconservative_terms, equations, + volume_integral::VolumeIntegralAdaptive{VolumeIntegralWeakForm, + VolumeIntegralFD, + Indicator}, + dg::DGSEM, + cache) where { + VolumeIntegralFD <: + VolumeIntegralFluxDifferencing, + Indicator <: AbstractIndicator} + @unpack volume_integral_default, volume_integral_stabilized, indicator = volume_integral + + # Calculate decision variable + decision = @trixi_timeit timer() "integral selector" indicator(u, mesh, equations, + dg, cache) + + @threaded for element in eachelement(dg, cache) + stabilized_version = decision[element] + + # TODO: Generalize/Dispatch or introduce yet sub-functions of the volume integrals + if stabilized_version + flux_differencing_kernel!(du, u, element, mesh, + nonconservative_terms, equations, + volume_integral_stabilized.volume_flux, + dg, cache) + else + weak_form_kernel!(du, u, element, mesh, + nonconservative_terms, equations, + dg, cache) + end + end + + return nothing +end + function calc_volume_integral!(du, u, mesh, have_nonconservative_terms, equations, volume_integral::VolumeIntegralPureLGLFiniteVolume, diff --git a/src/solvers/dgsem/dgsem.jl b/src/solvers/dgsem/dgsem.jl index 979f2d0577d..07aecc524c8 100644 --- a/src/solvers/dgsem/dgsem.jl +++ b/src/solvers/dgsem/dgsem.jl @@ -73,7 +73,8 @@ end Base.summary(io::IO, dg::DGSEM) = print(io, "DGSEM(polydeg=$(polydeg(dg)))") # `compute_u_mean` used in: -# (Stage-) Callbacks `EntropyBoundedLimiter` and `PositivityPreservingLimiterZhangShu` +# `IndicatorEntropyViolation` and the (stage-) limiters/callbacks +# `PositivityPreservingLimiterZhangShu` and `EntropyBoundedLimiter`. # positional arguments `mesh` and `cache` passed in to match signature of 2D/3D functions @inline function compute_u_mean(u::AbstractArray{<:Any, 3}, element, diff --git a/src/solvers/dgsem_tree/indicators.jl b/src/solvers/dgsem_tree/indicators.jl index 59847739c43..3a763228961 100644 --- a/src/solvers/dgsem_tree/indicators.jl +++ b/src/solvers/dgsem_tree/indicators.jl @@ -5,8 +5,6 @@ @muladd begin #! format: noindent -abstract type AbstractIndicator end - function create_cache(typ::Type{IndicatorType}, semi) where {IndicatorType <: AbstractIndicator} create_cache(typ, mesh_equations_solver_cache(semi)...) @@ -264,4 +262,122 @@ function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorMax) summary_box(io, "IndicatorMax", setup) end end + +@doc raw""" + IndicatorEntropyViolation(basis; + entropy_function=entropy, threshold=1e-9) + +This indicator checks the "entropy-violation" observed from timestep ``n \to n+1`` by +computing the difference in entropy ``\eta`` (determined by the `entropy_function`) as +```math +\Delta \eta_m^{(n+1)} = \eta^{(n+1)}(\bar{\boldsymbol{u}}_m) - \eta^{(n)}(\bar{\boldsymbol{u}}_m) +``` +for every cell `m` in the mesh. +The `entropy_function` is evaluated at the mean state ``\bar{\boldsymbol{u}}`` and defaults +to [`entropy`](@ref). + +The indicator returns for each cell ``m`` a boolean value indicating +whether the entropy growth/violation ``\Delta \eta_m^{(n+1)}`` exceeds the `threshold`. +Thus, for the mathematical entropy (default for `entropy`) which should not grow globally(!) +over the course of a simulation, this can be used to identify troubled cells. + +Supposed to be used in conjunction with [`VolumeIntegralAdaptive`](@ref) which then selects a +more advanced/ (entropy) stable volume integral for the troubled cells. + +!!! note + This indicator is not implemented as an AMR indicator, i.e., it is currently not + possible to employ this as the `indicator` in [`ControllerThreeLevel`](@ref). +""" +struct IndicatorEntropyViolation{EntropyFunction, RealT <: Real, Cache <: NamedTuple} <: + AbstractIndicator + entropy_function::EntropyFunction + threshold::RealT + cache::Cache +end + +# this method is used when the indicator is constructed as for adaptive volume integrals +function IndicatorEntropyViolation(basis; entropy_function = entropy, threshold = 1e-9) + cache = create_cache(IndicatorEntropyViolation, basis) + IndicatorEntropyViolation{typeof(entropy_function), + typeof(threshold), + typeof(cache)}(entropy_function, + threshold, cache) +end + +# this method is used when the indicator is constructed as for the adaptive volume integral +function create_cache(::Type{IndicatorEntropyViolation}, basis::LobattoLegendreBasis) + entropy_old = Vector{real(basis)}() + alpha = Vector{Bool}() + + return (; alpha, entropy_old) +end + +function Base.show(io::IO, indicator::IndicatorEntropyViolation) + @nospecialize indicator # reduce precompilation time + + print(io, "IndicatorEntropyViolation(") + print(io, "entropy_function=", indicator.entropy_function, ")") + print(io, ", threshold=", indicator.threshold, ")") +end + +function Base.show(io::IO, ::MIME"text/plain", indicator::IndicatorEntropyViolation) + @nospecialize indicator # reduce precompilation time + + if get(io, :compact, false) + show(io, indicator) + else + setup = [ + "entropy function" => indicator.entropy_function, + "threshold" => indicator.threshold + ] + summary_box(io, "IndicatorEntropyViolation", setup) + end +end + +function get_element_variables!(element_variables, indicator::IndicatorEntropyViolation, + ::VolumeIntegralAdaptive) + element_variables[:indicator_integral_stabilized] = Int.(indicator.cache.alpha) + return nothing +end + +# Dimension agnostic, i.e., valid for all 1D, 2D, and 3D meshes +function (indicator_entropy_violation::IndicatorEntropyViolation)(u, + mesh, equations, + dg::DGSEM, cache; + kwargs...) + @unpack alpha, entropy_old = indicator_entropy_violation.cache + resize!(alpha, nelements(dg, cache)) + @unpack entropy_function, threshold = indicator_entropy_violation + + # Beginning of simulation or after AMR: Need to compute `entropy_old` for every element + if length(entropy_old) != nelements(dg, cache) + resize!(entropy_old, nelements(dg, cache)) + + @threaded for element in eachelement(dg, cache) + u_mean = compute_u_mean(u, element, mesh, equations, dg, cache) + + # Compute entropy of the mean state + entropy_old[element] = entropy_function(u_mean, equations) + + alpha[element] = true # Be conservative: Use stabilized volume integral everywhere + end + else + @threaded for element in eachelement(dg, cache) + u_mean = compute_u_mean(u, element, mesh, equations, dg, cache) + + # Compute entropy of the mean state + entropy_element = entropy_function(u_mean, equations) + + # Check if entropy growth exceeds threshold + if entropy_element - entropy_old[element] > threshold + alpha[element] = true + else + alpha[element] = false + end + entropy_old[element] = entropy_element # Update `entropy_old` + end + end + + return alpha +end end # @muladd diff --git a/test/test_tree_1d_euler.jl b/test/test_tree_1d_euler.jl index 614dcc1b370..c844d34c6ac 100644 --- a/test/test_tree_1d_euler.jl +++ b/test/test_tree_1d_euler.jl @@ -340,6 +340,31 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_euler_adaptive_volume_integral.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_adaptive_volume_integral.jl"), + l2=[ + 0.11757168400166873, + 0.15258734107183258, + 0.43935902483141936 + ], + linf=[ + 0.19722297480582962, + 0.2575224502675185, + 0.7426620853498593 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) + + # Test/cover `show` + @test_nowarn show(stdout, indicator) + @test_nowarn show(IOContext(IOBuffer(), :compact => true), MIME"text/plain"(), + indicator) + @test_nowarn show(IOContext(IOBuffer(), :compact => false), MIME"text/plain"(), + volume_integral) +end + @trixi_testset "test_quasi_1D_entropy" begin using Trixi: CompressibleEulerEquationsQuasi1D, CompressibleEulerEquations1D, entropy, SVector diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl index da5305f0217..8a29db1a95d 100644 --- a/test/test_tree_2d_euler.jl +++ b/test/test_tree_2d_euler.jl @@ -837,6 +837,28 @@ end @test_allocations(Trixi.rhs!, semi, sol, 100) end +@trixi_testset "elixir_euler_warm_bubble_adaptive_integral.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_warm_bubble_adaptive_integral.jl"), + l2=[ + 0.00013825807605215672, + 0.020825709335780908, + 0.03327139463903476, + 31.40577277353277 + ], + linf=[ + 0.00168084149525638, + 0.15700603506542876, + 0.3344044463224682, + 330.7915698685101 + ], + tspan=(0.0, 10.0), + initial_refinement_level=4) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 100) +end + # Coverage test for all initial conditions @testset "Compressible Euler: Tests for initial conditions" begin @trixi_testset "elixir_euler_vortex.jl one step with initial_condition_constant" begin