diff --git a/examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl b/examples/dgmulti_2d/elixir_euler_kelvin_helmholtz_instability_adaptive_vol_int.jl new file mode 100644 index 00000000000..e3845ae78d5 --- /dev/null +++ b/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); diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl index e720d283ffa..21f0f308320 100644 --- a/src/callbacks_step/analysis_dgmulti.jl +++ b/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 diff --git a/src/solvers/dgmulti.jl b/src/solvers/dgmulti.jl deleted file mode 100644 index 363d91b5a4c..00000000000 --- a/src/solvers/dgmulti.jl +++ /dev/null @@ -1,17 +0,0 @@ -# includes solver files for DGMulti solvers -include("dgmulti/types.jl") -include("dgmulti/dg.jl") -include("dgmulti/flux_differencing_gauss_sbp.jl") -include("dgmulti/flux_differencing.jl") - -# integration of SummationByPartsOperators.jl -include("dgmulti/sbp.jl") - -# specialization of DGMulti to specific equations -include("dgmulti/flux_differencing_compressible_euler.jl") - -# shock capturing -include("dgmulti/shock_capturing.jl") - -# parabolic terms for DGMulti solvers -include("dgmulti/dg_parabolic.jl") diff --git a/src/solvers/dgmulti/dgmulti.jl b/src/solvers/dgmulti/dgmulti.jl new file mode 100644 index 00000000000..fabdf01e8ae --- /dev/null +++ b/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") diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl index 56d4e1f3752..a9b6813ccf3 100644 --- a/src/solvers/dgmulti/flux_differencing.jl +++ b/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, diff --git a/src/solvers/dgmulti/types.jl b/src/solvers/dgmulti/types.jl index 110432d254c..627ef809dfe 100644 --- a/src/solvers/dgmulti/types.jl +++ b/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, diff --git a/src/solvers/dgmulti/volume_integral_adaptive.jl b/src/solvers/dgmulti/volume_integral_adaptive.jl new file mode 100644 index 00000000000..0b29fdc3174 --- /dev/null +++ b/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 diff --git a/src/solvers/dgsem/special_volume_integrals.jl b/src/solvers/dgsem/special_volume_integrals.jl index 19ab8be7952..16cb0250599 100644 --- a/src/solvers/dgsem/special_volume_integrals.jl +++ b/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, diff --git a/src/solvers/solvers.jl b/src/solvers/solvers.jl index 8988a42bdab..5e242bb84bd 100644 --- a/src/solvers/solvers.jl +++ b/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 diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 4e00dd4859f..251b003522f 100644 --- a/test/test_dgmulti_2d.jl +++ b/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"),