creation of a vector containing the solution of the isentropic scheme and the...

creation of a vector containing the solution of the isentropic scheme and the explicit values of p and u

src/scheme/ImplicitAcousticSolver.cpp

@@ -42,10 +42,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
   BoundaryConditionList m_boundary_condition_list;
   const MeshType& m_mesh;
 
-  Vector<double> m_U;
   NodeValuePerCell<const Rdxd> m_Ajr;
   NodeValue<const bool> m_is_boundary_node;
 
+  CellValue<const Rd> m_predicted_u;
+  CellValue<const double> m_predicted_p;
+
   CellValue<const bool> m_is_implicit_cell;
   CellValue<const int> m_implicit_cell_index;
   size_t m_number_of_implicit_cells;
@@ -706,8 +708,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
                          const double& dt)
     : m_solver_type{solver_type},
       m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)},
-      m_mesh{*p_mesh},
-      m_U{0}
+      m_mesh{*p_mesh}
   {
     m_is_implicit_cell = [&] {
       CellValue<bool> is_implicit_cell(m_mesh.connectivity());
@@ -806,7 +807,25 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
 
     } while ((norm_inf_sol > 1e-14 * norm_inf_Un) and (nb_iter < 10000));
 
-    m_U = Uk;
+    m_predicted_u = [&] {
+      CellValue<Rd> predicted_u = copy(u.cellValues());
+      for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+        if (m_is_implicit_cell[cell_id]) {
+          predicted_u[cell_id] = Rd{Uk[mapU(0, cell_id)]};
+        }
+      }
+      return predicted_u;
+    }();
+
+    m_predicted_p = [&] {
+      CellValue<double> predicted_p = copy(p.cellValues());
+      for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+        if (m_is_implicit_cell[cell_id]) {
+          predicted_p[cell_id] = -Uk[mapP(cell_id)];
+        }
+      }
+      return predicted_p;
+    }();
   }
 
  public:
@@ -838,7 +857,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
           for (size_t j = 0; j < node_to_cell.size(); ++j) {
             const CellId J       = node_to_cell[j];
             const unsigned int R = node_local_number_in_its_cells[j];
-            br += m_Ajr(J, R) * Rd{m_U[mapU(0, J)]} - m_U[mapP(J)] * Cjr(J, R);
+            br += m_Ajr(J, R) * m_predicted_u[J] + m_predicted_p[J] * Cjr(J, R);
           }
 
           b[r] = br;
@@ -896,7 +915,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
           const auto& cell_nodes = cell_to_node_matrix[j];
 
           for (size_t r = 0; r < cell_nodes.size(); ++r) {
-            computed_Fjr(j, r) = m_Ajr(j, r) * (Rd{m_U[mapU(0, j)]} - ur[cell_nodes[r]]) - (m_U[mapP(j)] * Cjr(j, r));
+            computed_Fjr(j, r) = m_Ajr(j, r) * (m_predicted_u[j] - ur[cell_nodes[r]]) + (m_predicted_p[j] * Cjr(j, r));
           }
         });