diff --git a/src/solvers/dg.jl b/src/solvers/dg.jl
index 0e4d667fbc2..a0fb0d95079 100644
--- a/src/solvers/dg.jl
+++ b/src/solvers/dg.jl
@@ -12,6 +12,17 @@ function get_element_variables!(element_variables, u, mesh, equations,
     nothing
 end
 
+# Function to define "element variables" for the SaveSolutionCallback. It does
+# nothing by default, but can be specialized for certain mesh types. For instance,
+# parallel meshes output the mpi rank as an "element variable".
+function get_element_variables!(element_variables, mesh, dg, cache)
+    nothing
+end
+
+# Function to define "element variables" for the SaveSolutionCallback. It does
+# nothing by default, but can be specialized for certain volume integral types. 
+# For instance, shock capturing volume integrals output the blending factor
+# as an "element variable".
 function get_node_variables!(node_variables, mesh, equations,
                              volume_integral::AbstractVolumeIntegral, dg, cache)
     nothing
@@ -427,6 +438,7 @@ Base.summary(io::IO, dg::DG) = print(io, "DG(" * summary(dg.basis) * ")")
 function get_element_variables!(element_variables, u, mesh, equations, dg::DG, cache)
     get_element_variables!(element_variables, u, mesh, equations, dg.volume_integral,
                            dg, cache)
+    get_element_variables!(element_variables, mesh, dg, cache)
 end
 
 function get_node_variables!(node_variables, mesh, equations, dg::DG, cache)
diff --git a/src/solvers/dgsem_tree/dg_parallel.jl b/src/solvers/dgsem_tree/dg_parallel.jl
index c614fe0d0e6..3b33c9c33a5 100644
--- a/src/solvers/dgsem_tree/dg_parallel.jl
+++ b/src/solvers/dgsem_tree/dg_parallel.jl
@@ -5,6 +5,16 @@
 @muladd begin
 #! format: noindent
 
+# Function to output the mpi rank for visualization
+function get_element_variables!(element_variables,
+                                mesh::Union{ParallelTreeMesh, ParallelP4estMesh,
+                                            ParallelT8codeMesh},
+                                dg, cache)
+    element_variables[:mpi_rank] = ones(real(dg), nelements(mesh, dg, cache)) *
+                                   mpi_rank()
+    return nothing
+end
+
 # Initialize MPI data structures. This works for both the
 # `TreeMesh` and the `P4estMesh` and is dimension-agnostic.
 function init_mpi_data_structures(mpi_neighbor_interfaces, mpi_neighbor_mortars, n_dims,