Add implicit acoustic scheme structure

just renamed explicit acoustic scheme and plugged it

src/main.cpp

+
 #include <PugsUtils.hpp>
 
 #include <rang.hpp>
 
 #include <Connectivity.hpp>
 
-#include <AcousticSolver.hpp>
 #include <BoundaryCondition.hpp>
+#include <ImplicitAcousticSolver.hpp>
 #include <Mesh.hpp>
 
 #include <VTKWriter.hpp>
@@ -111,7 +112,7 @@ main(int argc, char* argv[])
 
       unknowns.initializeSod();
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      ImplicitAcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
       using Rd = TinyVector<MeshType::Dimension>;
 
@@ -231,7 +232,7 @@ main(int argc, char* argv[])
 
       unknowns.initializeSod();
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      ImplicitAcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
       const CellValue<const double>& Vj = mesh_data.Vj();
 
@@ -279,7 +280,7 @@ main(int argc, char* argv[])
       std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
                 << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+      method_cost_map["ImplicitAcousticSolverWithMesh"] = timer.seconds();
       break;
     }
     case 3: {
@@ -338,7 +339,7 @@ main(int argc, char* argv[])
 
       unknowns.initializeSod();
 
-      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      ImplicitAcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
       const CellValue<const double>& Vj = mesh_data.Vj();
 
@@ -385,7 +386,7 @@ main(int argc, char* argv[])
       std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
                 << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+      method_cost_map["ImplicitAcousticSolverWithMesh"] = timer.seconds();
       break;
     }
     }

src/scheme/ImplicitAcousticSolver.hpp

+#ifndef IMPLICIT_ACOUSTIC_SOLVER_HPP
+#define IMPLICIT_ACOUSTIC_SOLVER_HPP
+
+#include <rang.hpp>
+
+#include <ArrayUtils.hpp>
+
+#include <BlockPerfectGas.hpp>
+#include <PugsAssert.hpp>
+
+#include <BoundaryCondition.hpp>
+#include <FiniteVolumesEulerUnknowns.hpp>
+#include <Mesh.hpp>
+#include <MeshData.hpp>
+#include <SparseMatrixDescriptor.hpp>
+#include <TinyMatrix.hpp>
+#include <TinyVector.hpp>
+
+#include <ItemValueUtils.hpp>
+#include <Messenger.hpp>
+#include <SubItemValuePerItem.hpp>
+
+#include <iostream>
+
+template <typename MeshData>
+class ImplicitAcousticSolver
+{
+  using MeshType     = typename MeshData::MeshType;
+  using UnknownsType = FiniteVolumesEulerUnknowns<MeshData>;
+
+  MeshData& m_mesh_data;
+  const MeshType& m_mesh;
+  const typename MeshType::Connectivity& m_connectivity;
+  const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
+
+  constexpr static size_t Dimension = MeshType::Dimension;
+
+  using Rd  = TinyVector<Dimension>;
+  using Rdd = TinyMatrix<Dimension>;
+
+ private:
+  PUGS_INLINE
+  const CellValue<const double>
+  computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
+  {
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) { m_rhocj[j] = rhoj[j] * cj[j]; });
+    return m_rhocj;
+  }
+
+  PUGS_INLINE
+  const NodeValue<const double>
+  computeRhoCr(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
+  {
+    const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId& r) {
+        const auto& node_cells = node_to_cell_matrix[r];
+        const size_t nb_cells  = node_cells.size();
+        double rhocr           = 0;
+        for (size_t J = 0; J < nb_cells; ++J) {
+          CellId j = node_cells[J];
+          rhocr += rhoj[j] * cj[j];
+        }
+        m_rhocr[r] = rhocr * (1. / nb_cells);
+      });   // r[0] left border node, r[1] right border node
+    // then it fills in order
+    return m_rhocr;
+  }
+
+  PUGS_INLINE void
+  computeAjr(const CellValue<const double>& rhocj,
+             const NodeValuePerCell<const Rd>& Cjr,
+             const NodeValuePerCell<const double>& /* ljr */,
+             const NodeValuePerCell<const Rd>& njr)
+  {
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) {
+        const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
+        const double& rho_c    = rhocj[j];
+        for (size_t r = 0; r < nb_nodes; ++r) {
+          m_Ajr(j, r) = tensorProduct(rho_c * Cjr(j, r), njr(j, r));
+        }
+      });
+  }
+
+  PUGS_INLINE
+  const NodeValue<const Rdd>
+  computeAr(const NodeValuePerCell<const Rdd>& Ajr)
+  {
+    const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
+
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId& r) {
+        Rdd sum                                    = zero;
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
+
+        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];
+          sum += Ajr(J, R);
+        }
+        m_Ar[r] = sum;
+      });
+
+    return m_Ar;
+  }
+
+  PUGS_INLINE
+  const NodeValue<const Rd>
+  computeBr(const NodeValuePerCell<Rdd>& Ajr,
+            const NodeValuePerCell<const Rd>& Cjr,
+            const CellValue<const Rd>& uj,
+            const CellValue<const double>& pj)
+  {
+    const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
+
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId& r) {
+        Rd& br                                     = m_br[r];
+        br                                         = zero;
+        const auto& node_to_cell                   = node_to_cell_matrix[r];
+        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
+        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 += Ajr(J, R) * uj[J] + pj[J] * Cjr(J, R);
+        }
+      });
+
+    return m_br;
+  }
+
+  void
+  applyBoundaryConditions()
+  {
+    for (const auto& handler : m_boundary_condition_list) {
+      switch (handler.boundaryCondition().type()) {
+      case BoundaryCondition::normal_velocity: {
+        std::cerr << __FILE__ << ':' << __LINE__ << ": normal_velocity BC NIY\n";
+        std::exit(0);
+        break;
+      }
+      case BoundaryCondition::velocity: {
+        std::cerr << __FILE__ << ':' << __LINE__ << ": velocity BC NIY\n";
+        std::exit(0);
+        break;
+      }
+      case BoundaryCondition::pressure: {
+        // const PressureBoundaryCondition& pressure_bc
+        //   = dynamic_cast<const
+        //   PressureBoundaryCondition&>(handler.boundaryCondition());
+        std::cerr << __FILE__ << ':' << __LINE__ << ": pressure BC NIY\n";
+        std::exit(0);
+        break;
+      }
+      case BoundaryCondition::symmetry: {
+        const SymmetryBoundaryCondition<Dimension>& symmetry_bc =
+          dynamic_cast<const SymmetryBoundaryCondition<Dimension>&>(handler.boundaryCondition());
+        const Rd& n = symmetry_bc.outgoingNormal();
+
+        const Rdd I   = identity;
+        const Rdd nxn = tensorProduct(n, n);
+        const Rdd P   = I - nxn;
+
+        const Array<const NodeId>& node_list = symmetry_bc.nodeList();
+        parallel_for(
+          symmetry_bc.numberOfNodes(), PUGS_LAMBDA(const int& r_number) {
+            const NodeId r = node_list[r_number];
+
+            m_Ar[r] = P * m_Ar[r] * P + nxn;
+            m_br[r] = P * m_br[r];
+          });
+        break;
+      }
+      }
+    }
+  }
+
+  NodeValue<Rd>
+  computeUr(const NodeValue<const Rdd>& Ar, const NodeValue<const Rd>& br)
+  {
+    inverse(Ar, m_inv_Ar);
+    const NodeValue<const Rdd> invAr = m_inv_Ar;
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId& r) { m_ur[r] = invAr[r] * br[r]; });
+
+    return m_ur;
+  }
+
+  void
+  computeFjr(const NodeValuePerCell<Rdd>& Ajr,
+             const NodeValue<const Rd>& ur,
+             const NodeValuePerCell<const Rd>& Cjr,
+             const CellValue<const Rd>& uj,
+             const CellValue<const double>& pj)
+  {
+    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        for (size_t r = 0; r < cell_nodes.size(); ++r) {
+          m_Fjr(j, r) = Ajr(j, r) * (uj[j] - ur[cell_nodes[r]]) + pj[j] * Cjr(j, r);
+        }
+      });
+  }
+
+  void
+  inverse(const NodeValue<const Rdd>& A, NodeValue<Rdd>& inv_A) const
+  {
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId& r) { inv_A[r] = ::inverse(A[r]); });
+  }
+
+  PUGS_INLINE
+  void
+  computeImplicitFluxes(const CellValue<const double>& rhoj,
+                        const CellValue<const Rd>& uj,
+                        const CellValue<const double>& pj,
+                        const CellValue<const double>& cj,
+                        const NodeValuePerCell<const Rd>& Cjr,
+                        const NodeValuePerCell<const double>& ljr,
+                        const NodeValuePerCell<const Rd>& njr)
+  {
+    const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
+    computeAjr(rhocj, Cjr, ljr, njr);
+
+    NodeValuePerCell<const Rdd> Ajr = m_Ajr;
+    this->computeAr(Ajr);
+    this->computeBr(m_Ajr, Cjr, uj, pj);
+
+    this->applyBoundaryConditions();
+
+    synchronize(m_Ar);
+    synchronize(m_br);
+
+    NodeValue<Rd>& ur = m_ur;
+    ur                = computeUr(m_Ar, m_br);
+    computeFjr(m_Ajr, ur, Cjr, uj, pj);
+  }
+
+  NodeValue<Rd> m_br;
+  NodeValuePerCell<Rdd> m_Ajr;
+  NodeValue<Rdd> m_Ar;
+  NodeValue<Rdd> m_inv_Ar;
+  NodeValuePerCell<Rd> m_Fjr;
+  NodeValue<Rd> m_ur;
+  CellValue<double> m_rhocj;
+  CellValue<double> m_Vj_over_cj;
+  NodeValue<double> m_rhocr;
+
+ public:
+  ImplicitAcousticSolver(MeshData& mesh_data, const std::vector<BoundaryConditionHandler>& bc_list)
+    : m_mesh_data(mesh_data),
+      m_mesh(mesh_data.mesh()),
+      m_connectivity(m_mesh.connectivity()),
+      m_boundary_condition_list(bc_list),
+      m_br(m_connectivity),
+      m_Ajr(m_connectivity),
+      m_Ar(m_connectivity),
+      m_inv_Ar(m_connectivity),
+      m_Fjr(m_connectivity),
+      m_ur(m_connectivity),
+      m_rhocj(m_connectivity),
+      m_Vj_over_cj(m_connectivity),
+      m_rhocr(m_connectivity)
+  {
+    ;
+  }
+
+  PUGS_INLINE
+  double
+  acoustic_dt(const CellValue<const double>& Vj, const CellValue<const double>& cj) const
+  {
+    const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
+    const auto& cell_to_node_matrix           = m_mesh.connectivity().cellToNodeMatrix();
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double S = 0;
+        for (size_t r = 0; r < cell_nodes.size(); ++r) {
+          S += ljr(j, r);
+        }
+        m_Vj_over_cj[j] = 2 * Vj[j] / (S * cj[j]);
+      });
+
+    return min(m_Vj_over_cj);
+  }
+
+  void
+  computeNextStep(const double&, const double& dt, UnknownsType& unknowns)
+  {
+    CellValue<double>& rhoj = unknowns.rhoj();
+    CellValue<Rd>& uj       = unknowns.uj();
+    CellValue<double>& Ej   = unknowns.Ej();
+
+    CellValue<double>& ej = unknowns.ej();
+    CellValue<double>& pj = unknowns.pj();
+    CellValue<double>& cj = unknowns.cj();
+
+    const CellValue<const double>& Vj         = m_mesh_data.Vj();
+    const NodeValuePerCell<const Rd>& Cjr     = m_mesh_data.Cjr();
+    const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
+    const NodeValuePerCell<const Rd>& njr     = m_mesh_data.njr();
+
+    const NodeValue<const double> rhocr = computeRhoCr(rhoj, cj);
+
+    for (CellId j{0}; j < m_mesh.numberOfCells(); ++j) {
+      std::cout << "rho[" << j << "]=" << rhoj[j] << " c[" << j << "]=" << cj[j] << '\n';
+    }
+
+    for (NodeId r{0}; r < m_mesh.numberOfNodes(); ++r) {
+      std::cout << "rhocr[" << r << "]=" << rhocr[r] << '\n';
+    }
+
+    std::cout << "travailler\n";
+    std::exit(0);
+
+    computeImplicitFluxes(rhoj, uj, pj, cj, Cjr, ljr, njr);
+
+    const NodeValuePerCell<Rd>& Fjr = m_Fjr;
+    const NodeValue<const Rd> ur    = m_ur;
+    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+    const CellValue<const double> inv_mj = unknowns.invMj();
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd momentum_fluxes   = zero;
+        double energy_fluxes = 0;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          momentum_fluxes += Fjr(j, R);
+          energy_fluxes += (Fjr(j, R), ur[r]);
+        }
+        uj[j] -= (dt * inv_mj[j]) * momentum_fluxes;
+        Ej[j] -= (dt * inv_mj[j]) * energy_fluxes;
+      });
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) { ej[j] = Ej[j] - 0.5 * (uj[j], uj[j]); });
+
+    NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId& r) { mutable_xr[r] += dt * ur[r]; });
+    m_mesh_data.updateAllData();
+
+    const CellValue<const double> mj = unknowns.mj();
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId& j) { rhoj[j] = mj[j] / Vj[j]; });
+  }
+};
+
+#endif   // IMPLICIT_ACOUSTIC_SOLVER_HPP