`DGMulti` non-conservative terms appear to be broken

[jlchan]

I think the DGMulti treatment of all non-conservative terms is incorrect. In VolumeIntegralFluxDifferencing, I noticed that while conservative flux contributions are expected to be scaled by a factor of 2, non-conservative fluxes are not. This is not properly accounted for by DGMulti in hadamard_sum (which is used in the volume integral calculation).

The boundary terms incorporate the additional scaling of the non-conservative terms properly, but there may still be some issues related to the surface integral terms - if I fix the factor of 2 scaling issue in the volume integral, things still crash.

[jlchan]

This might explain some of the odd CI test failures with DGMulti for non-conservative equations we observed in previous PRs too.

[DanielDoehring]

Do you have found an example showcasing issues? This Alfven-wave test

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations
equations = IdealGlmMhdEquations2D(1.4)

surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell)
volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)

solver = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = GaussSBP(),
                 surface_integral = SurfaceIntegralWeakForm(surface_flux),
                 volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

initial_refinement_level = 3
N = 2^initial_refinement_level                 
cells_per_dimension = (N, N)
mesh = DGMultiMesh(solver, cells_per_dimension, periodicity = true, 
                   coordinates_min = (0.0, 0.0), coordinates_max = (sqrt(2), sqrt(2)))

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                     uEltype = real(solver))
alive_callback = AliveCallback(alive_interval = 10)

cfl = 0.8
stepsize_callback = StepsizeCallback(cfl = cfl)
glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
            dt = 1e-5, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep = false, callback = callbacks);

summary_callback() # print the timer summary

looks fine:

####################################################################################################
l2
rho                 rho_v1              rho_v2              rho_v3              rho_e               B1                  B2                  B3                  psi                 
error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       
1.87e-05  -         4.42e-06  -         4.42e-06  -         6.34e-06  -         7.21e-06  -         6.38e-06  -         6.38e-06  -         6.54e-06  -         4.40e-06  -         
5.28e-07  5.15      2.70e-07  4.03      2.70e-07  4.03      3.60e-07  4.14      4.46e-07  4.01      3.81e-07  4.06      3.81e-07  4.06      3.67e-07  4.16      4.80e-07  3.20      
2.72e-08  4.28      1.68e-08  4.01      1.68e-08  4.01      2.23e-08  4.01      2.98e-08  3.91      2.25e-08  4.08      2.25e-08  4.08      2.23e-08  4.04      4.41e-08  3.44      
1.64e-09  4.05      1.05e-09  4.00      1.05e-09  4.00      1.38e-09  4.02      2.04e-09  3.86      1.48e-09  3.92      1.48e-09  3.92      1.39e-09  4.01      1.03e-09  5.42      

mean      4.49      mean      4.01      mean      4.01      mean      4.06      mean      3.93      mean      4.02      mean      4.02      mean      4.07      mean      4.02      
----------------------------------------------------------------------------------------------------
linf
rho                 rho_v1              rho_v2              rho_v3              rho_e               B1                  B2                  B3                  psi                 
error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       error     EOC       
4.12e-05  -         1.07e-05  -         1.07e-05  -         1.46e-05  -         1.62e-05  -         1.36e-05  -         1.36e-05  -         1.44e-05  -         7.43e-06  -         
1.14e-06  5.17      6.17e-07  4.12      6.17e-07  4.12      7.85e-07  4.21      1.16e-06  3.80      8.10e-07  4.07      8.10e-07  4.07      8.08e-07  4.15      9.54e-07  2.96      
7.12e-08  4.00      4.03e-08  3.94      4.03e-08  3.94      4.92e-08  4.00      6.65e-08  4.13      4.71e-08  4.11      4.71e-08  4.11      4.93e-08  4.04      9.14e-08  3.38      
4.21e-09  4.08      2.54e-09  3.99      2.54e-09  3.99      3.11e-09  3.98      5.41e-09  3.62      3.30e-09  3.83      3.30e-09  3.83      3.10e-09  3.99      2.13e-09  5.42      

mean      4.42      mean      4.01      mean      4.01      mean      4.06      mean      3.85      mean      4.00      mean      4.00      mean      4.06      mean      3.92      
----------------------------------------------------------------------------------------------------

[jlchan]

I think this is because the non-conservative part of the ideal MHD equations is the divergence cleaning term. This can be incorrectly scaled without incurring a catastrophic error (it just amounts to a change in the cleaning speed).

If you try CompressibleEulerQuasi1D (see for example https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_1d_dgsem/elixir_euler_quasi_1d_source_terms.jl) with DGMulti, it fails while DGSEM + TreeMesh1D works.

[jlchan]

Nevermind; the scaling appears to be taken care of in https://github.com/jlchan/Trixi.jl/blob/a28a5600dfb2d41345287917742211ba274cf584/src/solvers/dgmulti/flux_differencing.jl#L487-L488 and https://github.com/jlchan/Trixi.jl/blob/a28a5600dfb2d41345287917742211ba274cf584/src/solvers/dgmulti/flux_differencing.jl#L544. I overlooked this when filing the issue.

Closing this issue.