Nonconservative boundary terms

[andrewwinters5000]

The split_form_kernel! found in solvers/dgsem_structured/dg_2d.jl uses the matrix derivative_split from dg.basis which is zero along the main diagonal. In essence, this "peels off" the two volume contributions at the [1,1] and [polydeg+1, polydeg+1] positions that would be present if one were to use derivative_matrix from dg.basis and moves them to behave more like boundary terms.

Currently for the nonconservative terms, the boundary flux must account for a contribution that depends only on interior values (which are in reality coming from the volume part) as well as part that couples the left and right states. See calc_boundary_flux! in solvers/dgsem_unstructured/dg_2d.jl for more information.

[andrewwinters5000]

The handling of the nonconservative terms in the volume / boundary will be the subject of a future PR. It will probably involve:

1. Writing an alternative split_form_kernel! that uses the standard derivative_matrix such that these diagonal (consistent) parts are no longer "peeled off" and they are directly computed as part of the volume.
2. Rederiving the nonconservative flux functions, e.g., flux_nonconservative_fjordholm_etal used for this surface / boundary coupling.
3. Some combination of the above two directions.

[ranocha]

If we used the standard derivative_matrix for nonconservative terms in the split_form_kernel!, we would also need to add additional terms at interfaces, wouldn't we? Thus, the way we handle this right now appears to be sane for SurfaceIntegralWeakForm. We could impose this explicitly in function signatures to clarify this aspect. Or is there a general alternative that works fine for both interfaces and boundaries?

[andrewwinters5000]

I need to look back at some experimenting I did on this subject for the shallow water equations. From what I remember, using the standard derivative_matrix would allow us to remove the "volume" contribution from flux_nonconservative_fjordholm_etal that only depends on u_ll and then the interfaces / boundaries for the nonconservative flux will only have the jump terms. For shallow water I remember this working, but I would need to test it out for the MHD equations.