Adaptive Volume Integral for DGMulti: WF & FD, Polynomial() & Affine

examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl

@@ -0,0 +1,86 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+volume_integral_weakform = VolumeIntegralWeakForm()
+volume_integral_fluxdiff = VolumeIntegralFluxDifferencing(flux_ranocha)
+
+# This indicator compares the entropy production of the weak form to the
+# true entropy evolution in that cell.
+# If the weak form does not increase entropy beyond `maximum_entropy_increase`,
+# we keep the weak form result. Otherwise, we switch to the stabilized/EC volume integral.
+indicator = IndicatorEntropyChange(maximum_entropy_increase = 5e-3)
+
+# Adaptive volume integral using the entropy production comparison indicator to perform the
+# stabilized/EC volume integral when needed.
+volume_integral = VolumeIntegralAdaptive(volume_integral_default = volume_integral_weakform,
+                                         volume_integral_stabilized = volume_integral_fluxdiff,
+                                         indicator = indicator)
+
+dg = DGMulti(polydeg = 3,
+             # `Tri()` and `Polynomial()` make flux differencing really(!) expensive
+             element_type = Tri(), approximation_type = Polynomial(),
+             surface_integral = SurfaceIntegralWeakForm(flux_hllc),
+             volume_integral = volume_integral)
+
+equations = CompressibleEulerEquations2D(1.4)
+
+"""
+    initial_condition_kelvin_helmholtz_instability(x, t, equations::CompressibleEulerEquations2D)
+
+A version of the classical Kelvin-Helmholtz instability based on
+- Andrés M. Rueda-Ramírez, Gregor J. Gassner (2021)
+  A Subcell Finite Volume Positivity-Preserving Limiter for DGSEM Discretizations
+  of the Euler Equations
+  [arXiv: 2102.06017](https://arxiv.org/abs/2102.06017)
+"""
+function initial_condition_kelvin_helmholtz_instability(x, t,
+                                                        equations::CompressibleEulerEquations2D)
+    # change discontinuity to tanh
+    # typical resolution 128^2, 256^2
+    # domain size is [-1,+1]^2
+    slope = 15
+    amplitude = 0.02
+    B = tanh(slope * x[2] + 7.5) - tanh(slope * x[2] - 7.5)
+    rho = 0.5 + 0.75 * B
+    v1 = 0.5 * (B - 1)
+    v2 = 0.1 * sin(2 * pi * x[1])
+    p = 1.0
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_kelvin_helmholtz_instability
+
+cells_per_dimension = (32, 32)
+mesh = DGMultiMesh(dg, cells_per_dimension; periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
+                                    boundary_conditions = boundary_condition_periodic)
+
+tspan = (0.0, 4.6) # stable time for limited entropy-increase adaptive volume integral
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+alive_callback = AliveCallback(alive_interval = 50)
+
+stepsize_callback = StepsizeCallback(cfl = 1.0)
+
+analysis_interval = 10
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval, uEltype = real(dg),
+                                     save_analysis = true,
+                                     analysis_errors = Symbol[],
+                                     extra_analysis_integrals = (entropy,))
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     solution_variables = cons2prim)
+
+callbacks = CallbackSet(summary_callback,
+                        alive_callback,
+                        stepsize_callback,
+                        analysis_callback,
+                        save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
+            dt = 1.0, ode_default_options()..., callback = callbacks);

src/callbacks_step/analysis_dgmulti.jl

@@ -134,10 +134,62 @@ function analyze(::Val{:linf_divb}, du, u, t,
     return linf_divB
 end
 
+# Calculate ∫_e (∂S/∂u ⋅ ∂u/∂t) dΩ_e where the result on element 'e' is kept in reference space
+# Note that ∂S/∂u = w(u) with entropy variables w.
+# This assumes that both du and u are already interpolated to the quadrature points
+function entropy_change_reference_element(du_values_local, u_values_local,
+                                          mesh::DGMultiMesh, equations,
+                                          dg::DGMulti, cache)
+    rd = dg.basis
+    @unpack Nq, wq = rd
+
+    # Compute entropy change for this element
+    dS_dt_elem = zero(eltype(first(du_values_local)))
+    for i in Base.OneTo(Nq) # Loop over quadrature points in the element
+        dS_dt_elem += dot(cons2entropy(u_values_local[i], equations),
+                          du_values_local[i]) * wq[i]
+    end
+
+    return dS_dt_elem
+end
+
+# calculate surface integral of func(u, normal_direction, equations) on the reference element.
+# For DGMulti, we loop over all faces of the element and integrate using face quadrature weights.
+# Restricted to `Polynomial` approximation type which requires interpolation to face quadrature nodes
+function surface_integral_reference_element(func::Func, u, element,
+                                            mesh::DGMultiMesh, equations,
+                                            dg::DGMulti,
+                                            cache, args...) where {Func}
+    rd = dg.basis
+    @unpack Nfq, wf, Vf = rd
+    md = mesh.md
+    @unpack nxyzJ = md
+
+    # Interpolate volume solution to face quadrature nodes for this element
+    @unpack u_face_local_threaded = cache
+    u_face_local = u_face_local_threaded[Threads.threadid()]
+    u_elem = view(u, :, element)
+    apply_to_each_field(mul_by!(Vf), u_face_local, u_elem)
+
+    surface_integral = zero(eltype(first(u)))
+    # Loop over all face nodes for this element
+    for i in 1:Nfq
+        # Get solution at this face node
+        u_node = u_face_local[i]
+
+        # Get face normal; nxyzJ stores components as (nxJ, nyJ, nxJ)
+        normal_direction = SVector(getindex.(nxyzJ, i, element))
+
+        # Multiply with face quadrature weight and accumulate
+        surface_integral += wf[i] * func(u_node, normal_direction, equations)
+    end
+
+    return surface_integral
+end
+
 function create_cache_analysis(analyzer, mesh::DGMultiMesh,
                                equations, dg::DGMulti, cache,
                                RealT, uEltype)
-    md = mesh.md
     return (;)
 end
 

src/solvers/dgmulti/dgmulti.jl

@@ -0,0 +1,22 @@
+# basic types and functions for DGMulti solvers
+include("types.jl")
+include("dg.jl")
+
+# flux differencing solver routines for DGMulti solvers
+include("flux_differencing_gauss_sbp.jl")
+include("flux_differencing.jl")
+
+# adaptive volume integral solver
+include("volume_integral_adaptive.jl")
+
+# integration of SummationByPartsOperators.jl
+include("sbp.jl")
+
+# specialization of DGMulti to specific equations
+include("flux_differencing_compressible_euler.jl")
+
+# shock capturing
+include("shock_capturing.jl")
+
+# parabolic terms for DGMulti solvers
+include("dg_parabolic.jl")

src/solvers/dgmulti/flux_differencing.jl

@@ -370,7 +370,7 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiff, RealT,
     # interpolate geometric terms to both quadrature and face values for curved meshes
     (; Vq, Vf) = dg.basis
     interpolated_geometric_terms = map(x -> [Vq; Vf] * x, mesh.md.rstxyzJ)
-    J = rd.Vq * md.J
+    J = Vq * md.J
 
     return (; md, Qrst_skew, VhP, Ph,
             invJ = inv.(J), dxidxhatj = interpolated_geometric_terms,

src/solvers/dgmulti/types.jl

@@ -22,9 +22,10 @@ const DGMultiWeakForm{ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxTyp
 const DGMultiFluxDiff{ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxType,
                                                       <:SurfaceIntegralWeakForm,
                                                       <:Union{VolumeIntegralFluxDifferencing,
-                                                              VolumeIntegralShockCapturingHGType}} where {
-                                                                                                          NDIMS
-                                                                                                          }
+                                                              VolumeIntegralShockCapturingHGType,
+                                                              VolumeIntegralAdaptiveEC_WF_DG}} where {
+                                                                                                      NDIMS
+                                                                                                      }
 
 const DGMultiFluxDiffSBP{ApproxType, ElemType} = DGMulti{NDIMS, ElemType, ApproxType,
                                                          <:SurfaceIntegralWeakForm,

src/solvers/dgmulti/volume_integral_adaptive.jl

@@ -0,0 +1,102 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+function create_cache(mesh::DGMultiMesh, equations,
+                      dg::DGMulti{NDIMS, ElemType, <:Polynomial,
+                                  <:SurfaceIntegralWeakForm,
+                                  <:VolumeIntegralAdaptive{<:IndicatorEntropyChange,
+                                                           <:VolumeIntegralWeakForm,
+                                                           <:VolumeIntegralFluxDifferencing}},
+                      RealT, uEltype) where {NDIMS, ElemType}
+    # Construct temporary solvers for each sub-integral to reuse the `create_cache` functions
+
+    # `VolumeIntegralAdaptive` for `DGMulti` currently limited to Weak Form & Flux Differencing combi
+    dg_WF = DG(dg.basis, dg.mortar, dg.surface_integral,
+               dg.volume_integral.volume_integral_default)
+    dg_FD = DG(dg.basis, dg.mortar, dg.surface_integral,
+               dg.volume_integral.volume_integral_stabilized)
+
+    cache_WF = create_cache(mesh, equations, dg_WF, RealT, uEltype)
+    cache_FD = create_cache(mesh, equations, dg_FD, RealT, uEltype)
+
+    # Set up structures required for `IndicatorEntropyChange`
+    rd = dg.basis
+    nvars = nvariables(equations)
+
+    # Required for entropy change computation (`entropy_change_reference_element`)
+    du_values = similar(cache_FD.u_values)
+
+    # Thread-local buffer for face interpolation, which is required
+    # for computation of entropy potential at interpolated face nodes
+    # (`surface_integral_reference_element`)
+    u_face_local_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nfq,), dg)
+                             for _ in 1:Threads.maxthreadid()]
+
+    return (; cache_FD...,
+            # Weak-form-specific fields for the default volume integral
+            weak_differentiation_matrices = cache_WF.weak_differentiation_matrices,
+            flux_threaded = cache_WF.flux_threaded,
+            rotated_flux_threaded = cache_WF.rotated_flux_threaded, # For non-affine meshes
+            # Required for `IndicatorEntropyChange`
+            du_values, u_face_local_threaded)
+end
+
+# version for affine meshes (currently only supported one for `VolumeIntegralAdaptive`)
+function calc_volume_integral!(du, u, mesh::DGMultiMesh,
+                               have_nonconservative_terms::False, equations,
+                               volume_integral::VolumeIntegralAdaptive{<:IndicatorEntropyChange},
+                               dg::DGMultiFluxDiff, cache)
+    @unpack volume_integral_default, volume_integral_stabilized = volume_integral
+    @unpack maximum_entropy_increase = volume_integral.indicator
+
+    # For weak form integral
+    @unpack u_values = cache
+
+    # For entropy production computation
+    rd = dg.basis
+    @unpack du_values = cache
+
+    # interpolate to quadrature points
+    apply_to_each_field(mul_by!(rd.Vq), u_values, u) # required for weak form trial
+
+    @threaded for e in eachelement(dg, cache)
+        # Try default volume integral first
+        volume_integral_kernel!(du, u, e, mesh,
+                                have_nonconservative_terms, equations,
+                                volume_integral_default, dg, cache)
+
+        # Interpolate `du` to quadrature points after WF integral for entropy production calculation
+        du_local = view(du, :, e)
+        du_values_local = view(du_values, :, e)
+        apply_to_each_field(mul_by!(rd.Vq), du_values_local, du_local) # required for entropy production calculation
+
+        # Compute entropy production of this volume integral
+        u_values_local = view(u_values, :, e)
+        dS_WF = -entropy_change_reference_element(du_values_local, u_values_local,
+                                                  mesh, equations,
+                                                  dg, cache)
+
+        dS_true = surface_integral_reference_element(entropy_potential, u, e,
+                                                     mesh, equations, dg, cache)
+
+        entropy_change = dS_WF - dS_true
+        if entropy_change > maximum_entropy_increase # Recompute using EC FD volume integral
+            # Reset default volume integral contribution.
+            # Note that this assumes that the volume terms are computed first,
+            # before any surface terms are added.
+            fill!(du_local, zero(eltype(du_local)))
+
+            # Recompute using stabilized volume integral. Note that the calculation of this volume integral requires the calculation of the entropy projection, which is done in `rhs!` specialized on the `DGMultiFluxDiff` solver type. 
+            volume_integral_kernel!(du, u, e, mesh,
+                                    have_nonconservative_terms, equations,
+                                    volume_integral_stabilized, dg, cache)
+        end
+    end
+
+    return nothing
+end
+end # @muladd

src/solvers/dgsem/special_volume_integrals.jl

@@ -7,6 +7,10 @@
 
 # This file contains some specialized volume integrals that require some indicators already to be defined.
 
+const VolumeIntegralAdaptiveEC_WF_DG = VolumeIntegralAdaptive{<:IndicatorEntropyChange,
+                                                              <:VolumeIntegralWeakForm,
+                                                              <:VolumeIntegralFluxDifferencing}
+
 """
     VolumeIntegralEntropyCorrection(indicator, 
                                     volume_integral_default, 

src/solvers/solvers.jl

@@ -20,5 +20,5 @@ end
 include("solvers_parabolic.jl")
 
 include("dg.jl")
-include("dgmulti.jl")
+include("dgmulti/dgmulti.jl")
 end # @muladd

test/test_dgmulti_2d.jl

@@ -354,6 +354,28 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl"),
+                        maximum_entropy_increase=0.0,
+                        tspan=(0.0, 0.2),
+                        l2=[
+                            0.05570371489805444,
+                            0.03299286402646503,
+                            0.05224508023471742,
+                            0.08011545946002244
+                        ],
+                        linf=[
+                            0.24323216643032874,
+                            0.1685158282708948,
+                            0.12357902305982191,
+                            0.26981068435988087
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_euler_rayleigh_taylor_instability.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_rayleigh_taylor_instability.jl"),