WIP debugging diamond scheme for heat equation [ci skip]

src/language/algorithms/HeatDiamondAlgorithm.cpp

@@ -34,7 +34,38 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
   std::shared_ptr m_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
   if constexpr (Dimension > 1) {
+    // if constexpr (Dimension == 2) {
+    //   NodeValue<TinyVector<Dimension>> new_coords{m_mesh->connectivity()};
+
+    //   new_coords[NodeId{0}] = TinyVector<2>{-1, -1};
+    //   new_coords[NodeId{1}] = TinyVector<2>{-1, 0};
+    //   new_coords[NodeId{2}] = TinyVector<2>{-1, 1};
+    //   new_coords[NodeId{3}] = TinyVector<2>{0, -1};
+    //   new_coords[NodeId{4}] = TinyVector<2>{0.2, 0.2};
+    //   new_coords[NodeId{5}] = TinyVector<2>{0, 1};
+    //   new_coords[NodeId{6}] = TinyVector<2>{1, -1};
+    //   new_coords[NodeId{7}] = TinyVector<2>{1, 0};
+    //   new_coords[NodeId{8}] = TinyVector<2>{1, 1};
+
+    //   m_mesh = std::make_shared<MeshType>(m_mesh->shared_connectivity(), new_coords);
+    // }
+
+    // for (NodeId node_id = 0; node_id < m_mesh->numberOfNodes(); ++node_id) {
+    //   std::cout << "x(" << node_id << ") = " << m_mesh->xr()[node_id] << '\n';
+    // }
+
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+    // {
+    //   const auto& cell_to_face_matrix = m_mesh->connectivity().cellToFaceMatrix();
+    //   CellId cell0                    = CellId{0};
+
+    //   const auto& ll = mesh_data.ll();
+
+    //   for (size_t i_face = 0; i_face < cell_to_face_matrix[cell0].size(); ++i_face) {
+    //     FaceId face_id = cell_to_face_matrix[cell0][i_face];
+    //     std::cout << "ll[" << face_id << "]=" << ll[face_id] << '\n';
+    //   }
+    // }
 
     CellValue<double> Tj =
       InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(T_id, mesh_data.xj());
@@ -188,9 +219,9 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
           for (size_t j_cell = 0; j_cell < primal_face_to_cell.size(); ++j_cell) {
             const CellId cell_id2 = primal_face_to_cell[j_cell];
             if (i_cell == j_cell) {
-              S(cell_id1, cell_id2) -= beta_l;
-            } else {
               S(cell_id1, cell_id2) += beta_l;
+            } else {
+              S(cell_id1, cell_id2) -= beta_l;
             }
           }
         }
@@ -209,19 +240,35 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
 
       NodeValue<const NodeId> dual_node_primal_node_id = [=]() {
         CellValue<NodeId> cell_ignored_id{m_mesh->connectivity()};
-        cell_ignored_id.fill(NodeId{std::numeric_limits<NodeId>::max()});
+        cell_ignored_id.fill(NodeId{std::numeric_limits<unsigned int>::max()});
+
         NodeValue<NodeId> node_primal_id{m_mesh->connectivity()};
 
         parallel_for(
-          m_mesh->numberOfCells(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
+          m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { node_primal_id[node_id] = node_id; });
 
         NodeValue<NodeId> computed_dual_node_primal_node_id{diamond_mesh->connectivity()};
 
         mapper->toDualNode(node_primal_id, cell_ignored_id, computed_dual_node_primal_node_id);
 
+        for (NodeId node_id = 0; node_id < diamond_mesh->numberOfNodes(); ++node_id) {
+          std::cout << node_id << " -> " << computed_dual_node_primal_node_id[node_id] << '\n';
+        }
+
         return computed_dual_node_primal_node_id;
       }();
 
+      // for (CellId cell_id = 0; cell_id < m_mesh->numberOfCells(); ++cell_id) {
+      //  std::cout << "xj[" << cell_id << "] = " << xj[cell_id] << "\n";
+      // }
+      // for (CellId dual_cell_id = 0; dual_cell_id < diamond_mesh->numberOfCells(); ++dual_cell_id) {
+      //  std::cout << "Vl[" << dual_cell_id << "] = " << dual_Vj[dual_cell_id] << "\n";
+      // }
+      // for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
+      //   std::cout << "alpha_l[" << face_id << "] = " << alpha_l[face_id] << "\n";
+      // }
+      // std::cout << "S(0,0) = " << S(0, 0) << '\n';
+
       const auto& dual_Cjr = diamond_mesh_data.Cjr();
 
       const auto& dual_cell_to_node_matrix = diamond_mesh->connectivity().cellToNodeMatrix();
@@ -231,19 +278,30 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
       const auto& face_local_numbers_in_their_cells = m_mesh->connectivity().faceLocalNumbersInTheirCells();
       const auto& primal_node_to_cell_matrix        = m_mesh->connectivity().nodeToCellMatrix();
 
+      const FaceValue<const TinyVector<Dimension>>& xl = mesh_data.xl();
+
       for (FaceId face_id = 0; face_id < m_mesh->numberOfFaces(); ++face_id) {
+        std::cout << "*** face_id -> " << face_id << "\n";
+        // std::cout << "    x_l -> " << xl[face_id] << "\n";
         if (face_to_cell_matrix[face_id].size() > 1) {
           const double alpha_face_id = alpha_l[face_id];
+          // std::cout << "    alpha_face_id -> " << alpha_face_id << "\n";
 
           for (size_t i_face_cell = 0; i_face_cell < face_to_cell_matrix[face_id].size(); ++i_face_cell) {
             CellId i_id = face_to_cell_matrix[face_id][i_face_cell];
-
+            // std::cout << "    cell_id_i_ -> " << i_id << "\n";
+            // std::cout << "    x_i -> " << xj[i_id] << "\n";
             const bool is_face_reversed_for_cell_i =
               primal_face_cell_is_reversed(i_id, face_local_numbers_in_their_cells(face_id, i_face_cell));
+            std::cout << "    i_is_reversed -> " << is_face_reversed_for_cell_i << "\n";
+
+            // for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
 
             for (size_t j_face_cell = 0; j_face_cell < face_to_cell_matrix[face_id].size(); ++j_face_cell) {
               CellId j_id         = face_to_cell_matrix[face_id][j_face_cell];
               CellId dual_cell_id = face_dual_cell_id[face_id];
+              // std::cout << "    cell_id_j -> " << j_id << "\n";
+              // std::cout << "    x_j -> " << xj[j_id] << "\n";
 
               for (size_t i_node = 0; i_node < primal_face_to_node_matrix[face_id].size(); ++i_node) {
                 const TinyVector<Dimension> nil = [&] {
@@ -253,9 +311,12 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
                     return primal_nlr(face_id, i_node);
                   }
                 }();
+                std::cout << "    nil -> " << nil << "[" << i_id << face_id << "]"
+                          << "\n";
 
                 NodeId face_node_id             = primal_face_to_node_matrix[face_id][i_node];
                 const auto& primal_node_to_cell = primal_node_to_cell_matrix[face_node_id];
+                std::cout << "    x_r -> " << xr[face_node_id] << "\n";
                 if (not primal_node_is_on_boundary[face_node_id]) {
                   const double weight_rj = [&] {
                     for (size_t i_cell = 0; i_cell < primal_node_to_cell.size(); ++i_cell) {
@@ -266,11 +327,14 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
                     Assert(false, "could not determine the weight");
                     return std::numeric_limits<double>::signaling_NaN();
                   }();
-                  for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
+                  for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size();
+                       ++i_dual_node) {
                     const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                    if (dual_node_id == face_node_id) {
+                    if (dual_node_primal_node_id[dual_node_id] == face_node_id) {
                       const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-                      S(i_id, j_id) += weight_rj * alpha_face_id * (nil, Clr);
+                      std::cout << "    Clr -> " << Clr << "\n";
+                      std::cout << "    w_rj -> " << weight_rj << "\n";
+                      S(i_id, j_id) -= weight_rj * alpha_face_id * (nil, Clr);
                     }
                   }
                 }   // else {
@@ -281,7 +345,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
           }
         }
       }
-
+      //      std::exit(0);
       CellValue<double> fj =
         InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id, mesh_data.xj());
 
@@ -299,36 +363,45 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
         b[cell_id]                      = fj[cell_id] * primal_Vj[cell_id];
         const auto& primal_cell_to_face = primal_cell_to_face_matrix[cell_id];
 
-        for (size_t face_id = 0; face_id < primal_cell_to_face.size(); ++face_id) {   // Num local
-          FaceId i_cell_face = primal_cell_to_face_matrix[cell_id][face_id];          // Num global
+        for (size_t i_cell_face = 0; i_cell_face < primal_cell_to_face.size(); ++i_cell_face) {
+          FaceId face_id = primal_cell_to_face_matrix[cell_id][i_cell_face];
 
           const bool is_face_reversed_for_cell_i = [=] {
-            for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {   // Num local
-              FaceId primal_face_id = primal_cell_to_face[i_face];                     // Num global
-              if (primal_face_id == i_cell_face) {
+            for (size_t i_face = 0; i_face < primal_cell_to_face.size(); ++i_face) {
+              FaceId primal_face_id = primal_cell_to_face[i_face];
+              if (primal_face_id == face_id) {
                 return primal_face_cell_is_reversed(cell_id, i_face);
               }
             }
+            Assert(false, "cannot find cell's face");
+            return false;
           }();
-          const auto& primal_face_to_node = primal_face_to_node_matrix[i_cell_face];
-          for (size_t i_node = 0; i_node < primal_face_to_node.size(); ++i_node) {   // Num local
+          const auto& primal_face_to_node = primal_face_to_node_matrix[face_id];
+          for (size_t i_node = 0; i_node < primal_face_to_node.size(); ++i_node) {
             const TinyVector<Dimension> nil = [=] {
               if (is_face_reversed_for_cell_i) {
-                return -primal_nlr(i_cell_face, i_node);
+                return -primal_nlr(face_id, i_node);
               } else {
-                return primal_nlr(i_cell_face, i_node);
+                return primal_nlr(face_id, i_node);
               }
             }();
 
-            CellId dual_cell_id = face_dual_cell_id[i_cell_face];   // Num global dual de la face
+            CellId dual_cell_id = face_dual_cell_id[face_id];
 
-            NodeId face_node_id = primal_face_to_node[i_node];   // Num global
+            NodeId face_node_id = primal_face_to_node[i_node];
             if (primal_node_is_on_boundary[face_node_id]) {
-              for (size_t i_dual_node = 0; i_dual_node < dual_Cjr.numberOfSubValues(dual_cell_id); ++i_dual_node) {
+              for (size_t i_dual_node = 0; i_dual_node < dual_cell_to_node_matrix[dual_cell_id].size(); ++i_dual_node) {
                 const NodeId dual_node_id = dual_cell_to_node_matrix[dual_cell_id][i_dual_node];
-                if (dual_node_id == face_node_id) {
+                if (dual_node_primal_node_id[dual_node_id] == face_node_id) {
                   const TinyVector<Dimension> Clr = dual_Cjr(dual_cell_id, i_dual_node);
-                  b[cell_id] -= alpha_l[i_cell_face] * gr[face_node_id] * (nil, Clr);
+                  b[cell_id] += alpha_l[face_id] * gr[face_node_id] * (nil, Clr);
+                  if (cell_id == 0) {
+                    std::cout << "***node_id" << face_node_id << "\n";
+                    std::cout << " *** alpha_l " << alpha_l[face_id] << "\n";
+                    std::cout << " *** Clr " << Clr << "\n";
+                    std::cout << " *** gr " << gr[face_node_id] << "\n";
+                    std::cout << " *** nil " << nil << "\n";
+                  }
                 }
               }
             }
@@ -337,8 +410,20 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
       }
 
       Vector<double> T{m_mesh->numberOfCells()};
-      T = 0;
-      BiCGStab{b, A, T, 1000, 1e-12};
+      T = 1;
+
+      for (size_t i = 0; i < b.size(); ++i) {
+        std::cout << "b[" << i << "]=" << b[i] << '\n';
+      }
+
+      Vector AT = A * T;
+
+      for (size_t i = 0; i < b.size(); ++i) {
+        std::cout << "AT[" << i << "]=" << AT[i] << '\n';
+      }
+
+      // PCG{b, A, A, T, 1000, 1e-12};
+      BiCGStab{b, A, T, 1000, 1e-15};
 
       CellValue<double> Temperature{m_mesh->connectivity()};
 
@@ -355,7 +440,7 @@ HeatDiamondScheme<Dimension>::HeatDiamondScheme(
       {
         VTKWriter vtk_writer("Temperature_" + std::to_string(Dimension), 0.01);
 
-        vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}}, 0,
+        vtk_writer.write(m_mesh, {NamedItemValue{"T", Temperature}, NamedItemValue{"Texact", Tj}}, 0,
                          true);   // forces last output
       }
     }