WIP: Additional analysis for multi-component Euler systems.

src/Trixi.jl

@@ -244,9 +244,9 @@ export cons2cons, cons2prim, prim2cons, cons2macroscopic, cons2state, cons2mean,
        cons2entropy, entropy2cons
 export density, pressure, density_pressure, velocity, global_mean_vars,
        equilibrium_distribution, waterheight, waterheight_pressure
-export entropy, energy_total, energy_kinetic, energy_internal,
-       energy_magnetic, cross_helicity, magnetic_field, divergence_cleaning_field,
-       enstrophy, vorticity
+export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic,
+       cross_helicity,
+       enstrophy, magnetic_field, divergence_cleaning_field, enstrophy_multi_euler, vorticity
 export lake_at_rest_error
 export ncomponents, eachcomponent
 

src/callbacks_step/analysis_dg2d.jl

@@ -246,10 +246,135 @@ function integrate(func::Func, u,
                                UnstructuredMesh2D, P4estMesh{2}, P4estMeshView{2},
                                T8codeMesh{2}},
                    equations, dg::DG, cache; normalize = true) where {Func}
-    integrate_via_indices(u, mesh, equations, dg, cache;
-                          normalize = normalize) do u, i, j, element, equations, dg
-        u_local = get_node_vars(u, equations, dg, i, j, element)
-        return func(u_local, equations)
+    @unpack boundaries = cache
+    m = methods(func)
+    if (m[1].nargs == 2) || (func == cons2cons)
+        return integrate_via_indices(u, mesh, equations, dg, cache;
+                                     normalize = normalize) do u, i, j, element,
+                                                               equations, dg
+            u_local = get_node_vars(u, equations, dg, i, j, element)
+
+            func(u_local, equations)
+        end
+    end
+    if (m[1].nargs == 4) && (func != cons2cons)
+        gradients_x = DGSpaceDerivative_WeakForm!(dg, cache, u, 1, equations, mesh)
+        gradients_y = DGSpaceDerivative_WeakForm!(dg, cache, u, 2, equations, mesh)
+        return integrate_via_indices(u, mesh, equations, dg, cache;
+                                     normalize = normalize) do u, i, j, element,
+                                                               equations, dg
+            u_local = get_node_vars(u, equations, dg, i, j, element)
+            gradients_local = Vector([gradients_x[:, i, j, element], gradients_y[:, i, j, element]])
+
+            func(u_local, gradients_local, equations)
+        end
+    end
+end
+
+# Return the derivatives of the interpolatied polynomial.
+function lagrange_derivatives(x, y)
+    n = length(x)
+    dydx = zeros(n)
+
+    for i in 1:n
+        for j in 1:n
+            if i != j
+                term = y[j] / (x[i] - x[j])
+                for k in 1:n
+                    if k != i && k != j
+                        term *= (x[i] - x[k]) / (x[j] - x[k])
+                    end
+                end
+                dydx[i] += term
+            end
+        end
+    end
+
+    return dydx
+end
+
+# Andrew's functions for computing the derivatives.
+function DGSpaceDerivative_WeakForm!(dg,
+                                     cache,
+                                     u,
+                                     direction::Int,
+                                     equations,
+                                     mesh)
+    # Get the required variables.
+    @unpack derivative_dhat, weights = dg.basis
+
+    # Compute the surface flux terms.
+    surface_flux_values = similar(u)
+    Trixi.calc_surface_integral!(surface_flux_values, u, mesh, equations, dg.surface_integral, dg, cache)
+
+    # Translations.
+    D = derivative_dhat
+    spA_Dhat = D
+
+    n_vars, Np, _, n_elements = size(u)
+    gradients = similar(u)
+
+    u_local = zeros((n_vars, Np, Np))
+
+    if direction == 1
+        for elem in 1:n_elements
+            # We work in primitive variables here.
+            for i in 1:Np
+                for j in 1:Np
+                    u_local[:, i, j] = cons2prim(u[:, i, j, elem], equations)
+                end
+            end
+            # Volume contributions (weak form)
+            for v in 1:n_vars
+                # Interpolate polynomial.
+                for j in 1:Np
+                    x = cache.elements.node_coordinates[1, :, j, elem]
+                    y = u_local[v, :, j]
+                    gradients[v, :, j, elem] .= lagrange_derivatives(x, y)
+                end
+
+#                 # Contract D in the ξ̂ (first spatial) direction
+#                 gradients[v, :, :, elem] .= D * u[v, :, :, elem]
+#                 # Average over element.
+
+                average = sum(gradients[v, :, :, elem]/Np^2)
+                gradients[v, :, :, elem] .= average
+            end
+        end
+    elseif direction == 2
+        for elem in 1:n_elements
+            # We work in primitive variables here.
+            for i in 1:Np
+                for j in 1:Np
+                    u_local[:, i, j] = cons2prim(u[:, i, j, elem], equations)
+                end
+            end
+            # Volume contributions (weak form)
+            for v in 1:n_vars
+                # Interpolate polynomial.
+                for i in 1:Np
+                    x = cache.elements.node_coordinates[2, i, :, elem]
+                    y = u_local[v, i, :]
+                    gradients[v, i, :, elem] .= lagrange_derivatives(x, y)
+                end
+
+#                 # Contract D in the ξ̂ (first spatial) direction
+#                 gradients[v, :, :, elem] .= u[v, :, :, elem] * D'
+
+                # Average over element.
+                average = sum(gradients[v, :, :, elem]/Np^2)
+                gradients[v, :, :, elem] .= average
+            end
+        end
+    end
+
+    return gradients
+end
+
+function compute_flux_array!(Flux, u, direction, equations)
+    Np, nvars = size(Flux)
+    for i in 1:Np
+        Flux[i, :] .= flux(u[i, :], 1, direction, equations)
     end
 end
 

src/equations/compressible_euler_multicomponent_2d.jl

@@ -979,4 +979,22 @@ end
     v = u[orientation] / rho
     return v
 end
+
+@inline function enstrophy_multi_euler(u, gradients,
+                           equations::CompressibleEulerMulticomponentEquations2D)
+    # Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v
+
+    omega = vorticity(u, gradients, equations)
+
+    return 0.5f0 * omega^2
+end
+
+@inline function vorticity(u, gradients,
+                           equations::CompressibleEulerMulticomponentEquations2D)
+    # Ensure that we have velocity `gradients` by way of the `convert_gradient_variables` function.
+    dv1dx, dv2dx = gradients[1][1:2]
+    dv1dy, dv2dy = gradients[2][1:2]
+
+    return dv2dx - dv1dy
+end
 end # @muladd