fix bugs and allow vectors and tensors

src/scheme/FluxingAdvectionSolver.cpp

@@ -126,10 +126,10 @@ FluxingAdvectionSolver<3>::computeFluxVolume() const
   throw NotImplementedError("ViensViensViens");
 }
 
-template <typename MeshType>
+template <typename MeshType, typename DataType>
 auto
 calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
-                     [[maybe_unused]] const FaceValue<double>& fluxing_volumes)
+                     [[maybe_unused]] const FaceValue<DataType>& fluxing_volumes)
 {
   constexpr size_t Dimension                               = MeshType::Dimension;
   const FaceValuePerCell<const bool> cell_face_is_reversed = old_mesh->connectivity().cellFaceIsReversed();
@@ -142,7 +142,7 @@ calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
       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)) {
+        if (!cell_face_is_reversed(cell_id, i_face)) {
           flux = -flux;
         }
         if (flux < 0) {
@@ -164,8 +164,9 @@ calculateRemapCycles(const std::shared_ptr<const MeshType>& old_mesh,
 
 template <typename MeshType, typename DataType>
 auto
-remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
-                  [[maybe_unused]] const FaceValue<double>& fluxing_volumes,
+remapUsingFluxing(const std::shared_ptr<const MeshType>& old_mesh,
+                  const std::shared_ptr<const MeshType>& new_mesh,
+                  const FaceValue<double>& fluxing_volumes,
                   const size_t num,
                   const DiscreteFunctionP0<MeshType::Dimension, const DataType>& old_q)
 {
@@ -173,17 +174,26 @@ remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
   //  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()));
+  DiscreteFunctionP0<Dimension, DataType> previous_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();
+  MeshData<Dimension>& old_mesh_data  = MeshDataManager::instance().getMeshData(*old_mesh);
+  const CellValue<const double> oldVj = old_mesh_data.Vj();
+  const CellValue<double> Vjstep(new_mesh->connectivity());
+  parallel_for(
+    new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+      Vjstep[cell_id] = oldVj[cell_id];
+      new_q[cell_id] *= oldVj[cell_id];
+    });
   for (size_t jstep = 0; jstep < num; ++jstep) {
-    std::cout << " step " << jstep << "\n";
+    // 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)) {
+          if (!cell_face_is_reversed(cell_id, i_face)) {
             flux = -flux;
           }
           const auto& face_to_cell = face_to_cell_matrix[face_id];
@@ -194,30 +204,49 @@ remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
           if (other_cell_id == cell_id) {
             other_cell_id = face_to_cell[1];
           }
-          DataType fluxed_q = old_q[cell_id];
+          DataType fluxed_q = previous_q[cell_id];
           if (flux > 0) {
-            fluxed_q = old_q[other_cell_id];
+            fluxed_q = previous_q[other_cell_id];
           }
+          Vjstep[cell_id] += flux / num;
           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";
       });
+    parallel_for(
+      new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        previous_q[cell_id] = 1 / Vjstep[cell_id] * new_q[cell_id];
+        //     std::cout << " old q " << old_q[cell_id] << " new q " << previous_q[cell_id] << "\n";
+        //     std::cout << " old vj " << oldVj[cell_id] << " new Vj " << Vjstep[cell_id] << "\n";
+      });
   }
+
+  MeshData<Dimension>& new_mesh_data  = MeshDataManager::instance().getMeshData(*new_mesh);
+  const CellValue<const double> newVj = new_mesh_data.Vj();
+  parallel_for(
+    new_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { new_q[cell_id] = 1 / newVj[cell_id] * new_q[cell_id]; });
+  // for (CellId cell_id = 0; cell_id < new_mesh->numberOfCells(); ++cell_id) {
+  //   if (abs(newVj[cell_id] - Vjstep[cell_id]) > 1e-15) {
+  //     std::cout << " cell " << cell_id << " newVj " << newVj[cell_id] << " Vjstep " << Vjstep[cell_id] << " diff "
+  //               << abs(newVj[cell_id] - Vjstep[cell_id]) << "\n";
+  //   }
+  // }
   return new_q;
 }
 
-template <typename MeshType>
+template <typename MeshType, typename DataType>
 auto
-remapUsingFluxing(const std::shared_ptr<const MeshType>& new_mesh,
+remapUsingFluxing([[maybe_unused]] const std::shared_ptr<const MeshType>& old_mesh,
+                  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)
+                  const DiscreteFunctionP0Vector<MeshType::Dimension, const DataType>& old_q)
 {
   constexpr size_t Dimension = MeshType::Dimension;
   //  const Connectivity<Dimension>& connectivity = new_mesh->connectivity();
 
-  DiscreteFunctionP0Vector<Dimension, double> new_q(new_mesh, copy(old_q.cellArrays()));
+  DiscreteFunctionP0Vector<Dimension, DataType> new_q(new_mesh, copy(old_q.cellArrays()));
 
   throw NotImplementedError("DiscreteFunctionP0Vector");
 
@@ -306,11 +335,11 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
 
 std::shared_ptr<const DiscreteFunctionVariant>
 FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
-                              const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variables)
+                              const std::shared_ptr<const DiscreteFunctionVariant>& remapped_variable)
 {
-  const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({remapped_variables});
+  const std::shared_ptr<const IMesh> old_mesh = getCommonMesh({remapped_variable});
 
-  if (not checkDiscretizationType({remapped_variables}, DiscreteFunctionType::P0)) {
+  if (not checkDiscretizationType({remapped_variable}, DiscreteFunctionType::P0)) {
     throw NormalError("acoustic solver expects P0 functions");
   }
 
@@ -330,9 +359,19 @@ FluxingAdvectionSolverHandler(const std::shared_ptr<const IMesh> new_mesh,
 
     FaceValue<double> fluxing_volumes = solver.computeFluxVolume();
     size_t number_of_cycles           = calculateRemapCycles(old_mesh0, fluxing_volumes);
-    return std::make_shared<DiscreteFunctionVariant>(
-      remapUsingFluxing(new_mesh0, fluxing_volumes, number_of_cycles,
-                        remapped_variables->get<DiscreteFunctionP0<Dimension, const double>>()));
+
+    DiscreteFunctionVariant new_variable = std::visit(
+      [&](auto&& variable) -> DiscreteFunctionVariant {
+        using DiscreteFunctionT = std::decay_t<decltype(variable)>;
+        if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
+          return remapUsingFluxing(old_mesh0, new_mesh0, fluxing_volumes, number_of_cycles, variable);
+        } else {
+          throw UnexpectedError("incompatible mesh types");
+        }
+      },
+      remapped_variable->discreteFunction());
+
+    return std::make_shared<DiscreteFunctionVariant>(new_variable);
   }
   case 3: {
     throw NotImplementedError("Fluxing advection solver not implemented in dimension 3");