Add MPI implementation for nonconservative equations on `TreeMesh` and `P4estMesh`

[JoshuaLampert]

When I run the elixir under examples/tree_2d_dgsem/elixir_shallowwater_source_terms.jl with MPI I get a MethodError: no method matching calc_mpi_interface_flux!. See the top of the stacktrace below (the whole stacktrace is just too long and imho the top is most interesting).

Stacktrace (on each rank)

```julia ERROR: LoadError: MethodError: no method matching calc_mpi_interface_flux!(::Array{Float64, 4}, ::TreeMesh{2, Trixi.ParallelTree{2}}, ::Static.True, ::ShallowWaterEquations2D{Float64}, ::SurfaceIntegralWeakForm{Tuple{FluxLaxFri[462/482] ypeof(max_abs_speed_naive)}, typeof(flux_nonconservative_fjordholm_etal)}}, ::DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 4, Matr ix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{Tuple{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}, typeof(flux_nonconservative_fjordholm_etal)}}, VolumeIntegralFluxDifferencing{Tuple{typeof(flux_wintermeyer_etal), typeof(flux_ nonconservative_wintermeyer_etal)}}}, ::NamedTuple{(:elements, :interfaces, :mpi_interfaces, :boundaries, :mortars, :mpi_mortars, :mpi_cache, :fstar_upper_threaded, :fstar_lower_threaded), Tuple{Trixi.ElementContainer2D{Float64, Float64 }, Trixi.InterfaceContainer2D{Float64}, Trixi.MPIInterfaceContainer2D{Float64}, Trixi.BoundaryContainer2D{Float64, Float64}, Trixi.L2MortarContainer2D{Float64}, Trixi.MPIL2MortarContainer2D{Float64}, Trixi.MPICache, Vector{StaticArraysC ore.MMatrix{4, 4, Float64, 16}}, Vector{StaticArraysCore.MMatrix{4, 4, Float64, 16}}}}) Closest candidates are: calc_mpi_interface_flux!(::Any, ::TreeMesh{2, <:Trixi.ParallelTree{2}}, ::Static.False, ::Any, ::Any, ::DG, ::Any) @ Trixi ~/work/julia/Trixi.jl/src/solvers/dgsem_tree/dg_2d_parallel.jl:682 calc_mpi_interface_flux!(::Any, ::P4estMesh{2, <:Real, <:Static.True}, ::Any, ::Any, ::Any, ::DG, ::Any) @ Trixi ~/work/julia/Trixi.jl/src/solvers/dgsem_p4est/dg_2d_parallel.jl:42 calc_mpi_interface_flux!(::Any, ::P4estMesh{3, <:Real, <:Static.True}, ::Any, ::Any, ::Any, ::DG, ::Any) @ Trixi ~/work/julia/Trixi.jl/src/solvers/dgsem_p4est/dg_3d_parallel.jl:141 ... Stacktrace: [1] macro expansion @ ~/work/julia/Trixi.jl/src/auxiliary/auxiliary.jl:272 [inlined] [2] rhs!(du::StrideArraysCore.PtrArray{Float64, 4, (1, 2, 3, 4), Tuple{Static.StaticInt{4}, Static.StaticInt{4}, Static.StaticInt{4}, Int64}, NTuple{4, Nothing}, NTuple{4, Static.StaticInt{1}}}, u::StrideArraysCore.PtrArray{Float64, 4 , (1, 2, 3, 4), Tuple{Static.StaticInt{4}, Static.StaticInt{4}, Static.StaticInt{4}, Int64}, NTuple{4, Nothing}, NTuple{4, Static.StaticInt{1}}}, t::Float64, mesh::TreeMesh{2, Trixi.ParallelTree{2}}, equations::ShallowWaterEquations2D{F loat64}, initial_condition::Function, boundary_conditions::Trixi.BoundaryConditionPeriodic, source_terms::Function, dg::DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trix i.LobattoLegendreMortarL2{Float64, 4, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{Tuple{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}, typeof(flux_nonconservative_fjordholm_etal)}}, VolumeIntegralFluxDifferencing{Tuple{t ypeof(flux_wintermeyer_etal), typeof(flux_nonconservative_wintermeyer_etal)}}}, cache::NamedTuple{(:elements, :interfaces, :mpi_interfaces, :boundaries, :mortars, :mpi_mortars, :mpi_cache, :fstar_upper_threaded, :fstar_lower_threaded), Tuple{Trixi.ElementContainer2D{Float64, Float64}, Trixi.InterfaceContainer2D{Float64}, Trixi.MPIInterfaceContainer2D{Float64}, Trixi.BoundaryContainer2D{Float64, Float64}, Trixi.L2MortarContainer2D{Float64}, Trixi.MPIL2MortarContainer2D {Float64}, Trixi.MPICache, Vector{StaticArraysCore.MMatrix{4, 4, Float64, 16}}, Vector{StaticArraysCore.MMatrix{4, 4, Float64, 16}}}}) @ Trixi ~/work/julia/Trixi.jl/src/solvers/dgsem_tree/dg_2d_parallel.jl:507 [3] macro expansion @ ~/work/julia/Trixi.jl/src/auxiliary/auxiliary.jl:272 [inlined] [434/482] [4] rhs!(du_ode::Vector{Float64}, u_ode::Vector{Float64}, semi::SemidiscretizationHyperbolic{TreeMesh{2, Trixi.ParallelTree{2}}, ShallowWaterEquations2D{Float64}, typeof(initial_condition_convergence_test), Trixi.BoundaryConditionPeri odic, typeof(source_terms_convergence_test), DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 4, Matrix{Float64}, Matrix{Float64}}, Su rfaceIntegralWeakForm{Tuple{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}, typeof(flux_nonconservative_fjordholm_etal)}}, VolumeIntegralFluxDifferencing{Tuple{typeof(flux_wintermeyer_etal), typeof(flux_nonconservative_wintermeyer_etal) }}}, NamedTuple{(:elements, :interfaces, :mpi_interfaces, :boundaries, :mortars, :mpi_mortars, :mpi_cache, :fstar_upper_threaded, :fstar_lower_threaded), Tuple{Trixi.ElementContainer2D{Float64, Float64}, Trixi.InterfaceContainer2D{Float 64}, Trixi.MPIInterfaceContainer2D{Float64}, Trixi.BoundaryContainer2D{Float64, Float64}, Trixi.L2MortarContainer2D{Float64}, Trixi.MPIL2MortarContainer2D{Float64}, Trixi.MPICache, Vector{StaticArraysCore.MMatrix{4, 4, Float64, 16}}, Ve ctor{StaticArraysCore.MMatrix{4, 4, Float64, 16}}}}}, t::Float64) @ Trixi ~/work/julia/Trixi.jl/src/semidiscretization/semidiscretization_hyperbolic.jl:297 [5] (::SciMLBase.Void{typeof(Trixi.rhs!)})(::Vector{Float64}, ::Vararg{Any}) @ SciMLBase ~/.julia/packages/SciMLBase/VdcHg/src/utils.jl:468 [6] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{typeof(Trixi.rhs!)}, arg1::Vector{Float64}, arg2::Vector{Float64}, arg3::SemidiscretizationHyperbolic{TreeMesh{2, Trixi.ParallelTree{2}}, ShallowWaterEquations2D{Float64} , typeof(initial_condition_convergence_test), Trixi.BoundaryConditionPeriodic, typeof(source_terms_convergence_test), DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi. LobattoLegendreMortarL2{Float64, 4, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{Tuple{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}, typeof(flux_nonconservative_fjordholm_etal)}}, VolumeIntegralFluxDifferencing{Tuple{typ eof(flux_wintermeyer_etal), typeof(flux_nonconservative_wintermeyer_etal)}}}, NamedTuple{(:elements, :interfaces, :mpi_interfaces, :boundaries, :mortars, :mpi_mortars, :mpi_cache, :fstar_upper_threaded, :fstar_lower_threaded), Tuple{Tri xi.ElementContainer2D{Float64, Float64}, Trixi.InterfaceContainer2D{Float64}, Trixi.MPIInterfaceContainer2D{Float64}, Trixi.BoundaryContainer2D{Float64, Float64}, Trixi.L2MortarContainer2D{Float64}, Trixi.MPIL2MortarContainer2D{Float64} , Trixi.MPICache, Vector{StaticArraysCore.MMatrix{4, 4, Float64, 16}}, Vector{StaticArraysCore.MMatrix{4, 4, Float64, 16}}}}}, arg4::Float64) @ FunctionWrappers ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:65 [7] macro expansion @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:137 [inlined] [8] do_ccall @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:125 [inlined] [9] FunctionWrapper @ ~/.julia/packages/FunctionWrappers/Q5cBx/src/FunctionWrappers.jl:144 [inlined] [10] _call @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:12 [inlined] [11] FunctionWrappersWrapper @ ~/.julia/packages/FunctionWrappersWrappers/9XR0m/src/FunctionWrappersWrappers.jl:10 [inlined] [12] ODEFunction @ ~/.julia/packages/SciMLBase/VdcHg/src/scimlfunctions.jl:2126 [inlined] [406/482] [13] macro expansion @ ~/.julia/packages/TimerOutputs/RsWnF/src/TimerOutput.jl:391 [inlined] [14] macro expansion @ ~/work/julia/Trixi.jl/src/callbacks_step/analysis.jl:292 [inlined] [15] macro expansion @ ~/work/julia/Trixi.jl/src/auxiliary/auxiliary.jl:272 [inlined] [16] (::AnalysisCallback{Trixi.LobattoLegendreAnalyzer{Float64, 7, SVector{7, Float64}, Matrix{Float64}}, Matrix{Float64}}, Tuple{typeof(Trixi.entropy_timederivative)}, SVector{4, Float64}, NamedTuple{(:u_local, :u_tmp1, :x_local, :x_tmp1), Tuple{Stride ArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{4}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 196}, Stri deArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{4}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 112}, St rideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{7}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 98}, S trideArraysCore.StaticStrideArray{Float64, 3, (1, 2, 3), Tuple{Static.StaticInt{2}, Static.StaticInt{7}, Static.StaticInt{4}}, Tuple{Nothing, Nothing, Nothing}, Tuple{Static.StaticInt{1}, Static.StaticInt{1}, Static.StaticInt{1}}, 56}}} }) [...] ```

[andrewwinters5000]

The problem is that for the MPI implementation on both TreeMesh and P4estMesh there only exists a calc_mpi_interface_flux! routine for the case when nonconservative_terms::False. The shallow water equations (as well as the MHD equations) need to have corresponding routines that computes the surface coupling that include the nonconservative terms. That is, there needs to be a version of this routine that would dispatch on nonconservative_terms::True and unpack the surface conservative and nonconservative flux functions out of the Tuple that is passed from the elixir.

[sloede]

Would you happen to know what we would need to add to support noncons terms for MPI? Is there anything fundamental that would have to be changed?

[andrewwinters5000]

Would you happen to know what we would need to add to support noncons terms for MPI? Is there anything fundamental that would have to be changed?

For TreeMesh it should be relatively straightforward as the nonconservative fluxes would simply need communicated. I believe the same is true for P4estMesh albeit a bit more complicated because of the nonuniqueness of the nonconservative terms and computing the mortar fluxes.

Addendum: I forgot that TreeMesh also has AMR 😬 . So the implementation procedure will likely be the same for both mesh types.

We might consider the IdealGLMMHDEquations to implement and test this new feature instead of ShallowWaterEquations. My thought process being that the shallow water equations with the current mortar implementation will not be well-balanced on any meshes with AMR (and thus should not be used).

I know that the new elixir from @JoshuaLampert in #1424 does not have AMR, but we need to be aware that users might misinterpret that the shallow water solver is ready to go on AMR meshes, which it is not.

[sloede]

In general, though, the ability to do non-conforming meshes and the ability to run MHD/SWE simulations in parallel are relatively orthogonal, aren't they? I am just checking to confirm that it could be a worthwhile endeavor to make the non-cons equations work in parallel on conforming meshes first, then add the ability to do mortars that remain well-balanced, and finally parallelize non-conforming meshes.

[andrewwinters5000]

In general, though, the ability to do non-conforming meshes and the ability to run MHD/SWE simulations in parallel are relatively orthogonal, aren't they?

Yes, absolutely! I just wanted to give a bit of a warning about the SWE and that testing against an MHD example might be better.

[JoshuaLampert]

Ok, thanks for the clarification @andrewwinters5000.

I am just checking to confirm that it could be a worthwhile endeavor to make the non-cons equations work in parallel on conforming meshes first, then add the ability to do mortars that remain well-balanced, and finally parallelize non-conforming meshes.

This sounds like a good plan.