Adaptive Volume Integral

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);

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);

src/Trixi.jl

@@ -262,6 +262,7 @@ export DG,
        VolumeIntegralFluxDifferencing,
        VolumeIntegralPureLGLFiniteVolume, VolumeIntegralPureLGLFiniteVolumeO2,
        VolumeIntegralShockCapturingHG, IndicatorHennemannGassner,
+       VolumeIntegralAdaptive, IndicatorEntropyViolation,
        VolumeIntegralUpwind,
        SurfaceIntegralWeakForm, SurfaceIntegralStrongForm,
        SurfaceIntegralUpwind,

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

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,

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,

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,