An optimized kernel for non-conservative fluxes

We've seen how the treatment of the geopotential term as non-conservative terms increases the robustness and the accuracy of the method. The same applies when someone wants to write down the vector invariant form as flux differencing. However the actual implementation is not optimal and could be greatly optimized by specifiying kernels when anti-symmetric fluxes are being used. I'm running the vector invariant form and I get a speed up of about 8 times when introducing these type of kernels: ```julia # Inlined function for interface flux computation for flux + nonconservative terms @inline function Trixi.calc_interface_flux!(surface_flux_values, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::True, equations::CompressibleEulerVectorInvariantEquations3D, surface_integral, dg::DG, cache, interface_index, normal_direction, primary_i_node_index, primary_j_node_index, primary_direction_index, primary_element_index, secondary_i_node_index, secondary_j_node_index, secondary_direction_index, secondary_element_index) @unpack u = cache.interfaces surface_flux, nonconservative_flux = surface_integral.surface_flux u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_i_node_index, primary_j_node_index, interface_index) flux_left, flux_right = surface_flux(u_ll, u_rr, normal_direction, equations) # Store the flux with nonconservative terms on the primary and secondary elements for v in eachvariable(equations) # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, primary_i_node_index, primary_j_node_index, primary_direction_index, primary_element_index] = flux_left[v] surface_flux_values[v, secondary_i_node_index, secondary_j_node_index, secondary_direction_index, secondary_element_index] = -flux_right[v] end return nothing end # inlined version of the boundary flux calculation along a physical interface @inline function Trixi.calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh::Union{P4estMesh{3}, T8codeMesh{3}}, nonconservative_terms::True, equations::CompressibleEulerVectorInvariantEquations3D, surface_integral, dg::DG, cache, i_index, j_index, k_index, i_node_index, j_node_index, direction_index, element_index, boundary_index) @unpack boundaries = cache @unpack node_coordinates, contravariant_vectors = cache.elements @unpack surface_flux = surface_integral # Extract solution data from boundary container u_inner = get_node_vars(boundaries.u, equations, dg, i_node_index, j_node_index, boundary_index) # Outward-pointing normal direction (not normalized) normal_direction = get_normal_direction(direction_index, contravariant_vectors, i_index, j_index, k_index, element_index) # Coordinates at boundary node x = get_node_coords(node_coordinates, equations, dg, i_index, j_index, k_index, element_index) # Call pointwise numerical flux functions for the conservative and nonconservative part # in the normal direction on the boundary flux = boundary_condition(u_inner, normal_direction, x, t, surface_flux, equations) # Copy flux to element storage in the correct orientation for v in eachvariable(equations) # Note the factor 0.5 necessary for the nonconservative fluxes based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs surface_flux_values[v, i_node_index, j_node_index, direction_index, element_index] = flux[v] end return nothing end ``` The same applies also for the volume flux, where we can exploit the same symmetric trick that's being used also for the non-conservative terms, where you just change a sign on the lower diagonal. These kernels just computes a single flux, which accounts for both conservative and non-conservative formulation. The mean values, especially if there are expensive means, they only get computed once and the anti-symmetry property is exploited by returning directly two fluxes, as a tuple, the left and the right. In case of pure conservative fluxes, these two fluxes, the right and left, coincide. An example of fluxes would be: ```julia @inline function flux_lmars_test(u_ll, u_rr, normal_direction::AbstractVector, equations::CompressibleEulerVectorInvariantEquations3D) #compute both conservative (f) and non conservative (g) fluxes. return SVector(f1, f2+ 0.5f0 * g2 * normal_direction[1], f3 + 0.5f0 * g3 * normal_direction[2], f4 + 0.5f0 * g4 * normal_direction[3], f5, 0), SVector(f1, f2 - 0.5f0 * g2 * normal_direction[1], f3 - 0.5f0 * g3 * normal_direction[2], f4 - 0.5f0 * g4 * normal_direction[3], f5, 0) end ``` This implementation is invasive in the current Trixi implementation, so the question is how we could introduce such a performance feature in the least invtrusive way?