Add remap step

src/scheme/FluxingAdvectionSolver.cpp

@@ -98,8 +98,10 @@ FluxingAdvectionSolver<2>::computeFluxVolume() const
   if (m_new_mesh->shared_connectivity() != m_old_mesh->shared_connectivity()) {
     throw NormalError("Old and new meshes must share the same connectivity");
   }
-  MeshDataType& old_mesh_data    = MeshDataManager::instance().getMeshData(*m_old_mesh);
-  MeshDataType& new_mesh_data    = MeshDataManager::instance().getMeshData(*m_new_mesh);
+  std::cout << " CARRE "
+            << "\n";
+  // MeshDataType& old_mesh_data    = MeshDataManager::instance().getMeshData(*m_old_mesh);
+  // MeshDataType& new_mesh_data    = MeshDataManager::instance().getMeshData(*m_new_mesh);
   const auto face_to_node_matrix = m_old_mesh->connectivity().faceToNodeMatrix();
   FaceValue<double> flux_volume(m_new_mesh->connectivity());
   parallel_for(
@@ -111,6 +113,8 @@ FluxingAdvectionSolver<2>::computeFluxVolume() const
       const Rd x3             = m_new_mesh->xr()[face_to_node[0]];
       TinyMatrix<2> M(x2[0] - x0[0], x3[0] - x1[0], x2[1] - x0[1], x3[1] - x1[1]);
       flux_volume[face_id] = 0.5 * det(M);
+      // std::cout << " x1 " << x1 << " x0 " << x0 << " x3 " << x3 << " x2 " << x2 << " flux volume "
+      //           << flux_volume[face_id] << "\n";
     });
   return flux_volume;
 }
@@ -122,17 +126,84 @@ FluxingAdvectionSolver<3>::computeFluxVolume() const
   throw NotImplementedError("ViensViensViens");
 }
 
+template <typename MeshType>
+auto
+calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
+                     [[maybe_unused]] const FaceValue<double>& fluxing_volumes)
+{
+  constexpr size_t Dimension                               = MeshType::Dimension;
+  const FaceValuePerCell<const bool> cell_face_is_reversed = old_mesh->connectivity().cellFaceIsReversed();
+  const auto cell_to_face_matrix                           = old_mesh->connectivity().cellToFaceMatrix();
+  const CellValue<double> total_negative_flux(old_mesh->connectivity());
+  total_negative_flux.fill(0);
+  parallel_for(
+    old_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      const auto& cell_to_face = cell_to_face_matrix[cell_id];
+      for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+        FaceId face_id = cell_to_face[i_face];
+        double flux    = fluxing_volumes[face_id];
+        if (cell_face_is_reversed(cell_id, i_face)) {
+          flux = -flux;
+        }
+        if (flux < 0) {
+          total_negative_flux[cell_id] += flux;
+        }
+      }
+      std::cout << " cell_id " << cell_id << " total_negative_flux " << total_negative_flux[cell_id] << "\n";
+    });
+  MeshData<Dimension>& mesh_data   = MeshDataManager::instance().getMeshData(*old_mesh);
+  const CellValue<const double> Vj = mesh_data.Vj();
+  const CellValue<size_t> ratio(old_mesh->connectivity());
+  parallel_for(
+    old_mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { ratio[cell_id] = std::ceil(abs(total_negative_flux[cell_id]) / Vj[cell_id]); });
+  size_t number_of_cycle = max(ratio);
+  std::cout << " number_of_cycle " << number_of_cycle << "\n";
+  return number_of_cycle;
+}
+
 template <typename MeshType, typename DataType>
 auto
 remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
                   [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
+                  const size_t num,
                   const DiscreteFunctionP0<MeshType::Dimension, const DataType>& old_q)
 {
   constexpr size_t Dimension = MeshType::Dimension;
   //  const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
-
+  const FaceValuePerCell<const bool> cell_face_is_reversed = new_mesh->connectivity().cellFaceIsReversed();
   DiscreteFunctionP0<Dimension, DataType> new_q(new_mesh, copy(old_q.cellValues()));
-
+  const auto cell_to_face_matrix = new_mesh->connectivity().cellToFaceMatrix();
+  const auto face_to_cell_matrix = new_mesh->connectivity().faceToCellMatrix();
+  for (size_t jstep = 0; jstep < num; ++jstep) {
+    std::cout << " step " << jstep << "\n";
+    parallel_for(
+      new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const auto& cell_to_face = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+          FaceId face_id = cell_to_face[i_face];
+          double flux    = fluxing_volumes[face_id];
+          if (cell_face_is_reversed(cell_id, i_face)) {
+            flux = -flux;
+          }
+          const auto& face_to_cell = face_to_cell_matrix[face_id];
+          if (face_to_cell.size() == 1) {
+            continue;
+          }
+          CellId other_cell_id = face_to_cell[0];
+          if (other_cell_id == cell_id) {
+            other_cell_id = face_to_cell[1];
+          }
+          DataType fluxed_q = old_q[cell_id];
+          if (flux > 0) {
+            fluxed_q = old_q[other_cell_id];
+          }
+          fluxed_q *= flux / num;
+          new_q[cell_id] += fluxed_q;
+        }
+        // std::cout << " old q " << old_q[cell_id] << " new q " << new_q[cell_id] << "\n";
+      });
+  }
   return new_q;
 }
 
@@ -140,6 +211,7 @@ template <typename MeshType>
 auto
 remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
                   [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
+                  [[maybe_unused]] const size_t num,
                   const DiscreteFunctionP0Vector<MeshType::Dimension, const double>& old_q)
 {
   constexpr size_t Dimension = MeshType::Dimension;
@@ -244,21 +316,7 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
 
   switch (old_mesh->dimension()) {
   case 1: {
-    constexpr size_t Dimension = 1;
-    using MeshType             = Mesh<Connectivity<Dimension>>;
-
-    const std::shared_ptr<const MeshType> old_mesh0 = std::dynamic_pointer_cast<const MeshType>(old_mesh);
-
-    const std::shared_ptr<const MeshType> new_mesh0 = std::dynamic_pointer_cast<const MeshType>(new_mesh);
-
-    // FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0);
-
-    FaceValue<double> fluxing_volumes(old_mesh0->connectivity());
-    fluxing_volumes.fill(0);
-
-    return std::make_shared<DiscreteFunctionVariant>(
-      remapUsingFluxing(new_mesh0, fluxing_volumes,
-                        remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>()));
+    throw NormalError("Not yet implemented in 1d");
   }
   case 2: {
     constexpr size_t Dimension = 2;
@@ -271,14 +329,10 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
     FluxingAdvectionSolver<Dimension> solver(old_mesh0, new_mesh0);
 
     FaceValue<double> fluxing_volumes = solver.computeFluxVolume();
-    //(new_mesh0->connectivity());
-    // fluxing_volumes.fill(0);
-
+    size_t number_of_cycles           = calculateRemapCycles(old_mesh0, fluxing_volumes);
     return std::make_shared<DiscreteFunctionVariant>(
-      remapUsingFluxing(new_mesh0, fluxing_volumes,
+      remapUsingFluxing(new_mesh0, fluxing_volumes, number_of_cycles,
                         remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>()));
-
-    // throw NotImplementedError("Fluxing advection solver not implemented in dimension 2");
   }
   case 3: {
     throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");