Fix mpi parallelism for flux advection solver

src/mesh/ItemArrayUtils.hpp

@@ -308,7 +308,7 @@ sum(const ItemArray<DataType, item_type, ConnectivityPtr>& item_value)
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
-synchronize(ItemArray<DataType, item_type, ConnectivityPtr>& item_array)
+synchronize(ItemArray<DataType, item_type, ConnectivityPtr> item_array)
 {
   static_assert(not std::is_const_v<DataType>, "cannot synchronize ItemArray of const data");
   if (parallel::size() > 1) {

src/mesh/ItemValueUtils.hpp

@@ -296,7 +296,7 @@ sum(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 void
-synchronize(ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
+synchronize(ItemValue<DataType, item_type, ConnectivityPtr> item_value)
 {
   static_assert(not std::is_const_v<DataType>, "cannot synchronize ItemValue of const data");
   if (parallel::size() > 1) {
@@ -309,7 +309,7 @@ synchronize(ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
 
 template <typename DataType, ItemType item_type, typename ConnectivityPtr>
 bool
-isSynchronized(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
+isSynchronized(ItemValue<DataType, item_type, ConnectivityPtr> item_value)
 {
   bool is_synchronized = true;
 

src/scheme/FluxingAdvectionSolver.cpp

@@ -3,6 +3,8 @@
 #include <language/utils/EvaluateAtPoints.hpp>
 #include <mesh/Connectivity.hpp>
 #include <mesh/IMesh.hpp>
+#include <mesh/ItemArrayUtils.hpp>
+#include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -150,16 +152,19 @@ template <>
 FaceValue<double>
 FluxingAdvectionSolver<1>::_computeAlgebraicFluxingVolume()
 {
-  Array<double> algebraic_fluxing_volumes{m_new_mesh->numberOfNodes()};
-  NodeValue<double> fluxing_volume(m_new_mesh->connectivity(), algebraic_fluxing_volumes);
+  Array<double> fluxing_volumes{m_new_mesh->numberOfNodes()};
+  NodeValue<double> nodal_fluxing_volume(m_new_mesh->connectivity(), fluxing_volumes);
   auto old_xr = m_old_mesh->xr();
   auto new_xr = m_new_mesh->xr();
 
   parallel_for(
     m_new_mesh->numberOfNodes(),
-    PUGS_LAMBDA(NodeId node_id) { fluxing_volume[node_id] = new_xr[node_id][0] - old_xr[node_id][0]; });
+    PUGS_LAMBDA(NodeId node_id) { nodal_fluxing_volume[node_id] = new_xr[node_id][0] - old_xr[node_id][0]; });
 
-  return FaceValue<double>{m_new_mesh->connectivity(), algebraic_fluxing_volumes};
+  FaceValue<double> algebraic_fluxing_volumes(m_new_mesh->connectivity(), fluxing_volumes);
+
+  synchronize(algebraic_fluxing_volumes);
+  return algebraic_fluxing_volumes;
 }
 
 template <>
@@ -185,6 +190,7 @@ FluxingAdvectionSolver<2>::_computeAlgebraicFluxingVolume()
       algebraic_fluxing_volume[face_id] = 0.5 * det(M);
     });
 
+  synchronize(algebraic_fluxing_volume);
   return algebraic_fluxing_volume;
 }
 
@@ -290,11 +296,12 @@ FluxingAdvectionSolver<Dimension>::_computeCycleNumber(FaceValue<double> fluxing
 
   MeshData<Dimension>& mesh_data   = MeshDataManager::instance().getMeshData(*m_old_mesh);
   const CellValue<const double> Vj = mesh_data.Vj();
-  const CellValue<size_t> ratio(m_old_mesh->connectivity());
+  CellValue<size_t> ratio(m_old_mesh->connectivity());
 
   parallel_for(
     m_old_mesh->numberOfCells(),
     PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(total_negative_flux[cell_id] / Vj[cell_id]); });
+  synchronize(ratio);
 
   size_t number_of_cycles = max(ratio);
 
@@ -379,6 +386,7 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
       }
     });
 
+  synchronize(new_q);
   old_q = new_q;
 }
 
@@ -411,6 +419,8 @@ FluxingAdvectionSolver<Dimension>::_remapAllQuantities()
         }
       });
 
+    synchronize(step_Vj);
+
     CellValue<double> inv_Vj(m_old_mesh->connectivity());
     parallel_for(
       m_new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { inv_Vj[cell_id] = 1 / step_Vj[cell_id]; });