WIP: ES volume integral

examples/2d/elixir_euler_kelvin_helmholtz_instability.jl

@@ -15,15 +15,17 @@ surface_flux = flux_lax_friedrichs
 volume_flux  = flux_chandrashekar
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
-indicator_sc = IndicatorHennemannGassner(equations, basis,
-                                         alpha_max=0.002,
-                                         alpha_min=0.0001,
-                                         alpha_smooth=true,
-                                         variable=density_pressure)
-volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
-                                                 volume_flux_dg=volume_flux,
-                                                 volume_flux_fv=surface_flux)
-solver = DGSEM(basis, surface_flux, volume_integral)
+#indicator_sc = IndicatorHennemannGassner(equations, basis,
+#                                         alpha_max=0.002,
+#                                         alpha_min=0.0001,
+#                                         alpha_smooth=true,
+#                                         variable=density_pressure)
+#volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+#                                                 volume_flux_dg=volume_flux,
+#                                                 volume_flux_fv=surface_flux)
+
+solver = DGSEM(3, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))
+#solver = DGSEM(basis, surface_flux, volume_integral)
 
 coordinates_min = (-0.5, -0.5)
 coordinates_max = ( 0.5,  0.5)

examples/2d/elixir_euler_source_terms.jl

@@ -10,7 +10,8 @@ equations = CompressibleEulerEquations2D(1.4)
 initial_condition = initial_condition_convergence_test
 
 surface_flux = flux_lax_friedrichs
-solver = DGSEM(3, surface_flux)
+volume_flux = flux_ranocha
+solver = DGSEM(3, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))
 
 coordinates_min = (0, 0)
 coordinates_max = (2, 2)

examples/3d/elixir_euler_taylor_green_vortex.jl

@@ -16,7 +16,7 @@ solver = DGSEM(3, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))
 coordinates_min = (-pi, -pi, -pi)
 coordinates_max = ( pi,  pi,  pi)
 mesh = TreeMesh(coordinates_min, coordinates_max,
-                initial_refinement_level=3,
+                initial_refinement_level=4,
                 n_cells_max=100_000)
 
 
@@ -31,8 +31,10 @@ ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
 
-analysis_interval = 100
-analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
+analysis_interval = 20
+analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
+                                     save_analysis=true,
+                                     extra_analysis_integrals=(Trixi.energy_kinetic,))
 
 alive_callback = AliveCallback(analysis_interval=analysis_interval)
 

src/solvers/dg/dg_2d.jl

@@ -24,6 +24,12 @@ function create_cache(mesh::TreeMesh{2}, equations::AbstractEquations{2},
   cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
   cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)
 
+  du_ec = Array{uEltype, 4}(undef, nvariables(equations), nnodes(dg), nnodes(dg),
+                            nelements(dg, cache))
+  du_cen = Array{uEltype, 4}(undef, nvariables(equations), nnodes(dg), nnodes(dg),
+                             nelements(dg, cache))
+  cache = (;cache..., du_ec, du_cen)
+
   return cache
 end
 
@@ -270,9 +276,45 @@ function calc_volume_integral!(du::AbstractArray{<:Any,4}, u,
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralFluxDifferencing,
                                dg::DGSEM, cache)
+
+  @unpack weights = dg.basis
+  @unpack du_ec, du_cen = cache
+
+  # du_ec = similar(du)
+  du_ec .= 0.0
+  # du_cen = similar(du)
+  du_cen .= 0.0
+
   @threaded for element in eachelement(dg, cache)
-    split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+    # compute volume integral with flux, and for comparison with central flux
+    split_form_kernel!(du_ec, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+    split_form_kernel!(du_cen, u, nonconservative_terms, equations, flux_central, dg, cache, element)
+
+    # compute entropy production of both volume integrals
+    delta_entropy = 0.0
+    for j in eachnode(dg), i in eachnode(dg)
+      du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element)
+      du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element)
+      w_node = cons2entropy(get_node_vars(u, equations, dg, i, j, element), equations)
+      delta_entropy +=weights[i]*weights[j]*dot(w_node,du_ec_node - du_cen_node)
+    end
+    if (delta_entropy < 0.0) 
+      for j in eachnode(dg), i in eachnode(dg)
+       du_cen_node = get_node_vars(du_cen, equations, dg, i, j, element)
+       add_to_node_vars!(du, du_cen_node, equations, dg, i, j, element)
+      end
+    else
+      for j in eachnode(dg), i in eachnode(dg)
+       du_ec_node = get_node_vars(du_ec, equations, dg, i, j, element)
+       add_to_node_vars!(du, du_ec_node, equations, dg, i, j, element)
+      end
+    end
   end
+
+  # Original code
+  # @threaded for element in eachelement(dg, cache)
+  #   split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+  # end
 end
 
 @inline function split_form_kernel!(du::AbstractArray{<:Any,4}, u,

src/solvers/dg/dg_3d.jl

@@ -24,6 +24,10 @@ function create_cache(mesh::TreeMesh{3}, equations::AbstractEquations{3},
   cache = (;cache..., create_cache(mesh, equations, dg.volume_integral, dg, uEltype)...)
   cache = (;cache..., create_cache(mesh, equations, dg.mortar, uEltype)...)
 
+  du_ec = Array{uEltype, 5}(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache))
+  du_cen = Array{uEltype, 5}(undef, nvariables(equations), nnodes(dg), nnodes(dg), nnodes(dg), nelements(dg, cache))
+  cache = (;cache..., du_ec, du_cen)
+
   return cache
 end
 
@@ -317,9 +321,43 @@ function calc_volume_integral!(du::AbstractArray{<:Any,5}, u,
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralFluxDifferencing,
                                dg::DGSEM, cache)
+  @unpack weights = dg.basis
+  @unpack du_ec, du_cen = cache
+
+  du_ec .= 0.0
+  du_cen .= 0.0
+
   @threaded for element in eachelement(dg, cache)
-    split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+    # compute volume integral with flux, and for comparison with central flux
+    split_form_kernel!(du_ec, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+    split_form_kernel!(du_cen, u, nonconservative_terms, equations, flux_central, dg, cache, element)
+    # split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+
+    # compute entropy production of both volume integrals
+    delta_entropy = 0.0
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+      du_ec_node = get_node_vars(du_ec, equations, dg, i, j, k, element)
+      du_cen_node = get_node_vars(du_cen, equations, dg, i, j, k, element)
+      w_node = cons2entropy(get_node_vars(u, equations, dg, i, j, k, element), equations)
+      delta_entropy +=weights[i]*weights[j]*weights[k]*dot(w_node,du_ec_node - du_cen_node)
+    end
+    if (delta_entropy < 0.0) 
+      for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+       du_cen_node = get_node_vars(du_cen, equations, dg, i, j, k, element)
+       add_to_node_vars!(du, du_cen_node, equations, dg, i, j, k, element)
+      end
+    else
+      for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+       du_ec_node = get_node_vars(du_ec, equations, dg, i, j, k, element)
+       add_to_node_vars!(du, du_ec_node, equations, dg, i, j, k, element)
+      end
+    end
   end
+
+  # Original code
+  # @threaded for element in eachelement(dg, cache)
+  #   split_form_kernel!(du, u, nonconservative_terms, equations, volume_integral.volume_flux, dg, cache, element)
+  # end
 end
 
 @inline function split_form_kernel!(du::AbstractArray{<:Any,5}, u,