Added first working version of MPI with non-conservative systems for …

…CONFORMING meshes (analysis routines need work)

examples/p4est_3d_dgsem/elixir_mhd_alfven_wave_conforming.jl

@@ -0,0 +1,64 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations3D(5 / 3)
+
+initial_condition = initial_condition_convergence_test
+
+volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
+solver = DGSEM(polydeg = 3,
+               surface_flux = (flux_hlle,
+                               flux_nonconservative_powell),
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0, -1.0, -1.0)
+coordinates_max = (1.0, 1.0, 1.0)
+
+# Create P4estMesh with 2 x 2 x 2 trees
+trees_per_dimension = (2, 2, 2)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 3, initial_refinement_level = 2,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 10
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+cfl = 1.0
+stepsize_callback = StepsizeCallback(cfl = cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback,
+                        glm_speed_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # 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

src/equations/ideal_glm_mhd_3d.jl

@@ -41,7 +41,7 @@ function varnames(::typeof(cons2prim), ::IdealGlmMhdEquations3D)
     ("rho", "v1", "v2", "v3", "p", "B1", "B2", "B3", "psi")
 end
 function default_analysis_integrals(::IdealGlmMhdEquations3D)
-    (entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
+    (entropy_timederivative,)#Val(:l2_divb), Val(:linf_divb)) # Temporarily deactivated because parallel analysis is not working
 end
 
 # Set initial conditions at physical location `x` for time `t`

src/solvers/dgsem_p4est/dg_3d_parallel.jl

@@ -267,6 +267,41 @@ end
     end
 end
 
+# Inlined version of the interface flux computation for non-conservative equations
+@inline function calc_mpi_interface_flux!(surface_flux_values,
+                                          mesh::Union{ParallelP4estMesh{3},
+                                                      ParallelT8codeMesh{3}},
+                                          nonconservative_terms::True, equations,
+                                          surface_integral, dg::DG, cache,
+                                          interface_index, normal_direction,
+                                          interface_i_node_index,
+                                          interface_j_node_index, local_side,
+                                          surface_i_node_index, surface_j_node_index,
+                                          local_direction_index, local_element_index)
+    @unpack u = cache.mpi_interfaces
+    surface_flux, nonconservative_flux = surface_integral.surface_flux
+
+    u_ll, u_rr = get_surface_node_vars(u, equations, dg,
+                                       interface_i_node_index, interface_j_node_index,
+                                       interface_index)
+
+    # Compute flux and non-conservative term for this side of the interface
+    if local_side == 1
+        flux_ = surface_flux(u_ll, u_rr, normal_direction, equations)
+        noncons_flux_ = 0.5 *
+                        nonconservative_flux(u_ll, u_rr, normal_direction, equations)
+    else # local_side == 2
+        flux_ = -surface_flux(u_ll, u_rr, -normal_direction, equations)
+        noncons_flux_ = -0.5 *
+                        nonconservative_flux(u_rr, u_ll, -normal_direction, equations)
+    end
+
+    for v in eachvariable(equations)
+        surface_flux_values[v, surface_i_node_index, surface_j_node_index,
+        local_direction_index, local_element_index] = flux_[v] + noncons_flux_[v]
+    end
+end
+
 function prolong2mpimortars!(cache, u,
                              mesh::Union{ParallelP4estMesh{3}, ParallelT8codeMesh{3}},
                              equations,