Add an operator *= to SparseMatrixDescriptor. Provide a parabolic heat solver

src/algebra/SparseMatrixDescriptor.hpp

@@ -34,16 +34,14 @@ class SparseMatrixDescriptor
       return m_id_value_map.size();
     }
 
-    const DataType&
-    operator[](const IndexType& j) const
+    const DataType& operator[](const IndexType& j) const
     {
       auto i_value = m_id_value_map.find(j);
       Assert(i_value != m_id_value_map.end());
       return i_value->second;
     }
 
-    DataType&
-    operator[](const IndexType& j)
+    DataType& operator[](const IndexType& j)
     {
       auto i_value = m_id_value_map.find(j);
       if (i_value != m_id_value_map.end()) {
@@ -100,6 +98,18 @@ class SparseMatrixDescriptor
     return m_row_array[i][j];
   }
 
+  SparseMatrixDescriptor&
+  operator*=(const double& lambda)
+  {
+    for (IndexType i_row = 0; i_row < m_row_array.size(); ++i_row) {
+      SparseRowDescriptor& row = m_row_array[i_row];
+      for (auto& [id, value] : row.m_id_value_map) {
+        value *= lambda;
+      }
+    }
+    return *this;
+  }
+
   const DataType&
   operator()(const IndexType& i, const IndexType& j) const
   {

src/language/algorithms/ParabolicHeat.cpp

@@ -31,7 +31,7 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
   const FunctionSymbolId& T_id,
   const FunctionSymbolId& kappa_id,
   const FunctionSymbolId& f_id,
-  const double& final_time,
+  const double& Tf,
   const double& dt)
 {
   using ConnectivityType = Connectivity<Dimension>;
@@ -252,7 +252,15 @@ ParabolicHeatScheme<Dimension>::ParabolicHeatScheme(
                         NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
                        0, true);   // forces last output
     }
-    ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face);
+    double time = 0;
+    Assert(dt > 0, "The time-step have to be positive!");
+    do {
+      double deltat = std::min(dt, Tf - time);
+      std::cout << "Current time = " << time << " time-step = " << deltat << " final time = " << Tf << "\n";
+      ScalarDiamondScheme<Dimension>(i_mesh, bc_descriptor_list, kappa_id, f_id, Temperature, Temperature_face, Tf,
+                                     deltat);
+      time += deltat;
+    } while (time < Tf && std::abs(time - Tf) > 1e-15);
     {
       Vector<double> error{mesh->numberOfCells()};
       CellValue<double> cell_error{mesh->connectivity()};

src/scheme/ScalarDiamondScheme.hpp

@@ -150,7 +150,9 @@ class ScalarDiamondScheme
                       const FunctionSymbolId& kappa_id,
                       const FunctionSymbolId& f_id,
                       CellValue<double>& Temperature,
-                      FaceValue<double>& Temperature_face)
+                      FaceValue<double>& Temperature_face,
+                      const double& Tf,
+                      const double& dt)
   {
     using ConnectivityType = Connectivity<Dimension>;
     using MeshType         = Mesh<ConnectivityType>;
@@ -583,6 +585,11 @@ class ScalarDiamondScheme
           return computed_alpha_l;
         }();
 
+        double lambda      = (Tf == 0) ? 0 : 1;
+        double time_factor = (lambda == 0) ? 1 : dt;
+        std::cout << " lambda = " << lambda << " time_factor =" << time_factor << "\n";
+        const CellValue<const double> primal_Vj = mesh_data.Vj();
+
         SparseMatrixDescriptor S{number_of_dof};
         for (FaceId face_id = 0; face_id < mesh->numberOfFaces(); ++face_id) {
           const auto& primal_face_to_cell = face_to_cell_matrix[face_id];
@@ -593,12 +600,13 @@ class ScalarDiamondScheme
             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_dof_number[cell_id1], cell_dof_number[cell_id2]) += beta_l;
+                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) +=
+                  (time_factor * beta_l + lambda * primal_Vj[cell_id1]);
                 if (primal_face_is_neumann[face_id]) {
                   S(face_dof_number[face_id], cell_dof_number[cell_id2]) -= beta_l;
                 }
               } else {
-                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= beta_l;
+                S(cell_dof_number[cell_id1], cell_dof_number[cell_id2]) -= time_factor * beta_l;
               }
             }
           }
@@ -636,7 +644,7 @@ class ScalarDiamondScheme
 
                   for (size_t j_cell = 0; j_cell < primal_node_to_cell_matrix[node_id].size(); ++j_cell) {
                     CellId j_id = primal_node_to_cell_matrix[node_id][j_cell];
-                    S(cell_dof_number[i_id], cell_dof_number[j_id]) -= w_rj(node_id, j_cell) * a_ir;
+                    S(cell_dof_number[i_id], cell_dof_number[j_id]) -= time_factor * w_rj(node_id, j_cell) * a_ir;
                     if (primal_face_is_neumann[face_id]) {
                       S(face_dof_number[face_id], cell_dof_number[j_id]) += w_rj(node_id, j_cell) * a_ir;
                     }
@@ -645,7 +653,7 @@ class ScalarDiamondScheme
                     for (size_t l_face = 0; l_face < node_to_face_matrix[node_id].size(); ++l_face) {
                       FaceId l_id = node_to_face_matrix[node_id][l_face];
                       if (primal_face_is_on_boundary[l_id]) {
-                        S(cell_dof_number[i_id], face_dof_number[l_id]) -= w_rl(node_id, l_face) * a_ir;
+                        S(cell_dof_number[i_id], face_dof_number[l_id]) -= time_factor * w_rl(node_id, l_face) * a_ir;
                         if (primal_face_is_neumann[face_id]) {
                           S(face_dof_number[face_id], face_dof_number[l_id]) += w_rl(node_id, l_face) * a_ir;
                         }
@@ -667,12 +675,12 @@ class ScalarDiamondScheme
           InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(f_id,
                                                                                                     mesh_data.xj());
 
-        const CellValue<const double> primal_Vj = mesh_data.Vj();
         CRSMatrix A{S};
         Vector<double> b{number_of_dof};
 
         for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
-          b[cell_dof_number[cell_id]] = fj[cell_id] * primal_Vj[cell_id];
+          b[cell_dof_number[cell_id]] =
+            (time_factor * fj[cell_id] + lambda * Temperature[cell_id]) * primal_Vj[cell_id];
         }
         // Dirichlet on b^L_D
         {
@@ -703,7 +711,8 @@ class ScalarDiamondScheme
                 const auto& value_list = bc.valueList();
                 for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
                   FaceId face_id = face_list[i_face];
-                  b[face_dof_number[face_id]] += mes_l[face_id] * value_list[i_face];   // sign
+                  b[face_dof_number[face_id]] +=
+                    mes_l[face_id] * value_list[i_face];   // + lambda * Temperature_face[face_id]);   // sign
                 }
               }
             },
@@ -735,7 +744,9 @@ template ScalarDiamondScheme<1>::ScalarDiamondScheme(
   const FunctionSymbolId&,
   const FunctionSymbolId&,
   CellValue<double>&,
-  FaceValue<double>&);
+  FaceValue<double>&,
+  const double&,
+  const double&);
 
 template ScalarDiamondScheme<2>::ScalarDiamondScheme(
   std::shared_ptr<const IMesh>,
@@ -743,7 +754,9 @@ template ScalarDiamondScheme<2>::ScalarDiamondScheme(
   const FunctionSymbolId&,
   const FunctionSymbolId&,
   CellValue<double>&,
-  FaceValue<double>&);
+  FaceValue<double>&,
+  const double&,
+  const double&);
 
 template ScalarDiamondScheme<3>::ScalarDiamondScheme(
   std::shared_ptr<const IMesh>,
@@ -751,4 +764,6 @@ template ScalarDiamondScheme<3>::ScalarDiamondScheme(
   const FunctionSymbolId&,
   const FunctionSymbolId&,
   CellValue<double>&,
-  FaceValue<double>&);
+  FaceValue<double>&,
+  const double&,
+  const double&);

tests/test_SparseMatrixDescriptor.cpp

@@ -196,4 +196,31 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
 )";
     REQUIRE(output.str() == expected_output);
   }
+  SECTION("Multiplication by scalar")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 2) = 11;
+    S(1, 1) = 1;
+    S(0, 2) += 2;
+    S(3, 3) = 5;
+    S(3, 1) = -3;
+    S(4, 1) = 1;
+    S(2, 2) = 4;
+    S(4, 4) = -2;
+    S *= 3;
+    const auto value_array = S.valueArray();
+
+    REQUIRE(value_array.size() == 9);
+
+    REQUIRE(value_array[0] == 0);
+    REQUIRE(value_array[1] == 15);
+    REQUIRE(value_array[2] == 3);
+    REQUIRE(value_array[3] == 33);
+    REQUIRE(value_array[4] == 12);
+    REQUIRE(value_array[5] == -9);
+    REQUIRE(value_array[6] == 15);
+    REQUIRE(value_array[7] == 3);
+    REQUIRE(value_array[8] == -6);
+  }
 }