Add code to test a 1d Lagrangian solver based on Glace formalism

src/main.cpp

@@ -8,6 +8,7 @@
 #include <Mesh.hpp>
 #include <BoundaryCondition.hpp>
 #include <AcousticSolver.hpp>
+#include <VoronoiSolver.hpp>
 
 #include <VTKWriter.hpp>
 
@@ -41,6 +42,128 @@ int main(int argc, char *argv[])
 
       std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
 
+      const bool use_voronoi
+          = [&] () {
+              char* voronoi_env = getenv("VORONOI");
+              return voronoi_env != nullptr;
+            } ();
+
+      if (use_voronoi) {
+        if (p_mesh->meshDimension() != 1) {
+          std::cerr << "higher dimension than 1 is not implemented for Voronoi\n";
+          std::abort();
+        }
+
+        {
+          std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
+          std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+          for (const auto& sym_boundary_name : sym_boundary_name_list){
+            std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+                = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+            SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+                = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+            bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+          }
+
+          using ConnectivityType = Connectivity1D;
+          using MeshType = Mesh<ConnectivityType>;
+          using MeshDataType = MeshData<MeshType>;
+          using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+          const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+          Timer timer;
+          timer.reset();
+          MeshDataType mesh_data(mesh);
+
+          std::vector<BoundaryConditionHandler> bc_list;
+          {
+            for (const auto& bc_descriptor : bc_descriptor_list) {
+              switch (bc_descriptor->type()) {
+                case BoundaryConditionDescriptor::Type::symmetry: {
+                  const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
+                      = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+                  for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList();
+                       ++i_ref_node_list) {
+                    const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list);
+                    const RefId& ref = ref_node_list.refId();
+                    if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                      SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+                          = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_node_list));
+                      std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+                      bc_list.push_back(BoundaryConditionHandler(bc));
+                    }
+                  }
+                  break;
+                }
+                default: {
+                  perr() << "Unknown BCDescription\n";
+                  std::exit(1);
+                }
+              }
+            }
+          }
+
+          UnknownsType unknowns(mesh_data);
+
+          unknowns.initializeSod();
+
+          VoronoiSolver<MeshDataType> voronoi_solver(mesh_data, unknowns, bc_list);
+
+          using Rd = TinyVector<MeshType::dimension>;
+
+          const CellValue<const double>& Vj = mesh_data.Vj();
+
+          const double tmax=0.2;
+          double t=0;
+
+          int itermax=std::numeric_limits<int>::max();
+          int iteration=0;
+
+          CellValue<double>& rhoj = unknowns.rhoj();
+          CellValue<double>& ej = unknowns.ej();
+          CellValue<double>& pj = unknowns.pj();
+          CellValue<double>& gammaj = unknowns.gammaj();
+          CellValue<double>& cj = unknowns.cj();
+
+          BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+
+          VTKWriter vtk_writer("mesh", 0.01);
+
+          while((t<tmax) and (iteration<itermax)) {
+            std::cout << "time " << t << ' ';
+            vtk_writer.write(mesh, t);
+            double dt = 0.4*voronoi_solver.voronoi_dt(Vj, cj);
+            if (t+dt>tmax) {
+              dt=tmax-t;
+            }
+            voronoi_solver.computeNextStep(t,dt, unknowns);
+            std::cout << "dt=" << dt << '\n';
+            block_eos.updatePandCFromRhoE();
+
+            t += dt;
+            ++iteration;
+          }
+          vtk_writer.write(mesh, t, true); // forces last output
+
+          pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
+                 << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+          method_cost_map["VoronoiSolverWithMesh"] = timer.seconds();
+
+          { // gnuplot output for density
+            const CellValue<const Rd>& xj = mesh_data.xj();
+            const CellValue<const double>& rhoj = unknowns.rhoj();
+            std::ofstream fout("rho");
+            for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+              fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+            }
+          }
+
+
+        }
+      } else {
         switch (p_mesh->meshDimension()) {
           case 1: {
             std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
@@ -343,7 +466,7 @@ int main(int argc, char *argv[])
             break;
           }
         }
-
+      }
       pout() << "* "  << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
 
     } else {

src/mesh/MeshData.hpp

@@ -179,7 +179,7 @@ class MeshData
               const auto& face_nodes = face_to_node_matrix[l];
 
 #warning should this lambda be replaced by a precomputed correspondance?
-              std::function local_node_number_in_cell
+              auto local_node_number_in_cell
                   = [&](const NodeId& node_number) {
                       for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
                         if (node_number == cell_nodes[i_node]) {

src/scheme/VoronoiSolver.hpp

+#ifndef VORONOI_SOLVER_HPP
+#define VORONOI_SOLVER_HPP
+
+#include <rang.hpp>
+
+#include <ArrayUtils.hpp>
+
+#include <BlockPerfectGas.hpp>
+#include <PastisAssert.hpp>
+
+#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+#include <Mesh.hpp>
+#include <MeshData.hpp>
+#include <FiniteVolumesEulerUnknowns.hpp>
+#include <BoundaryCondition.hpp>
+
+#include <SubItemValuePerItem.hpp>
+
+template<typename MeshData>
+class VoronoiSolver
+{
+  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:
+  PASTIS_INLINE
+  const CellValue<const double>
+  computeRhoCj(const CellValue<const double>& rhoj,
+               const CellValue<const double>& cj)
+  {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        m_rhocj[j] = rhoj[j]*cj[j];
+      });
+    return m_rhocj;
+  }
+
+  PASTIS_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(), PASTIS_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));
+        }
+      });
+  }
+
+  PASTIS_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(), PASTIS_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;
+  }
+
+  PASTIS_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(), PASTIS_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: {
+          perr() << __FILE__ << ':' << __LINE__  << ": normal_velocity BC NIY\n";
+          std::exit(0);
+          break;
+        }
+        case BoundaryCondition::velocity: {
+          perr() << __FILE__ << ':' << __LINE__  << ": velocity BC NIY\n";
+          std::exit(0);
+          break;
+        }
+        case BoundaryCondition::pressure: {
+          // const PressureBoundaryCondition& pressure_bc
+          //   = dynamic_cast<const PressureBoundaryCondition&>(handler.boundaryCondition());
+          perr() << __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(), PASTIS_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(), PASTIS_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(), PASTIS_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(), PASTIS_LAMBDA(const NodeId& r) {
+        inv_A[r] = ::inverse(A[r]);
+      });
+  }
+
+  PASTIS_INLINE
+  void computeExplicitFluxes(const NodeValue<const Rd>& xr,
+                             const CellValue<const Rd>& xj,
+                             const CellValue<const double>& rhoj,
+                             const CellValue<const Rd>& uj,
+                             const CellValue<const double>& pj,
+                             const CellValue<const double>& cj,
+                             const CellValue<const double>& Vj,
+                             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;
+    const NodeValue<const Rdd> Ar = computeAr(Ajr);
+    const NodeValue<const Rd> br = computeBr(m_Ajr, Cjr, uj, pj);
+
+    this->applyBoundaryConditions();
+
+    NodeValue<Rd>& ur = m_ur;
+    ur = computeUr(Ar, 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;
+
+  CellValue<Rd> m_xg;
+
+ public:
+  VoronoiSolver(MeshData& mesh_data,
+                 UnknownsType& unknowns,
+                 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_xg(m_connectivity)
+  {
+    ;
+  }
+
+  PASTIS_INLINE
+  double voronoi_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(), PASTIS_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 ReduceMin(m_Vj_over_cj);
+  }
+
+  void computeNextStep(const double& t, 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 Rd>& xj = m_mesh_data.xj();
+
+    static bool init_xg=false;
+    if (not init_xg) {
+      parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+          m_xg[j] = xj[j];
+        });
+      init_xg = true;
+    }
+
+    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 Rd>& xr = m_mesh.xr();
+
+    // computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr);
+
+    CellValue<Rd> uj_star(m_mesh.connectivity());
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        if ((j != 0) and ( j != m_mesh.numberOfCells()-1)) {
+          CellId jp1 = j+1;
+          CellId jm1 = j-1;
+          double rhoc_jp1 = rhoj[jp1]*cj[jp1];
+          double rhoc_jm1 = rhoj[jm1]*cj[jm1];
+          uj_star[j]
+              = (1./ (rhoc_jp1+rhoc_jm1)) * (rhoc_jp1 * uj[jp1] + rhoc_jm1*uj[jm1])
+              + (pj[jm1] - pj[jp1]) / (rhoc_jp1+rhoc_jm1);
+        } else {
+          uj_star[j] = zero;
+        }
+      });
+
+    // for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
+    //   std::cout << "uj_star[" << j << "]= " << uj_star[j] << "\n";
+    // }
+
+
+    CellValue<double> pj_star(m_mesh.connectivity());
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        if ((j > 0) and ( j < m_mesh.numberOfCells()-1)) {
+          CellId jm1 = j-1;
+          double rhoc_jm1 = rhoj[jm1]*cj[jm1];
+          pj_star[j] = rhoc_jm1*(uj[jm1]-uj_star[j])[0] + pj[jm1];
+        } else {
+          pj_star[j] = pj[j];
+        }
+      });
+
+    const CellValue<const double> inv_mj = unknowns.invMj();
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        if ((j > 0) and ( j < m_mesh.numberOfCells()-1)) {
+          CellId jp1 = j+1;
+          CellId jm1 = j-1;
+          uj[j] -= (dt*inv_mj[j]) * 0.5*(pj_star[jp1]-pj_star[jm1]);
+          Ej[j] -= (dt*inv_mj[j]) * 0.5*(pj_star[jp1]*uj_star[jp1][0]-pj_star[jm1]*uj_star[jm1][0]);
+        }
+      });
+
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        m_xg[j] += dt*uj_star[j];
+      });
+
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
+      });
+
+    const auto& node_to_cell_matrix
+        = m_connectivity.nodeToCellMatrix();
+
+    NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r){
+        const auto& node_to_cell = node_to_cell_matrix[r];
+
+        if (node_to_cell.size() == 2) {
+          Rd new_xr = zero;
+          for (size_t j=0; j<node_to_cell.size(); ++j) {
+            const CellId J = node_to_cell[j];
+            new_xr += m_xg[J];
+          }
+          mutable_xr[r] = 0.5*new_xr;
+        }
+      });
+    m_mesh_data.updateAllData();
+
+    const CellValue<const double> mj = unknowns.mj();
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+        rhoj[j] = mj[j]/Vj[j];
+      });
+
+
+    // for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
+    //   std::cout << "rho[" << j << "]= " << rhoj[j] << "\t ";
+    //   std::cout << "u[" << j << "]= " << uj[j] << "\t ";
+    //   std::cout << "E[" << j << "]= " << Ej[j] << "\t ";
+    //   std::cout << "e[" << j << "]= " << ej[j] << "\n";
+    // }
+    // std::cout << std::flush;
+    // std::cout << "enter value to continue: ";
+    // int i;
+    // std::cin >> i;
+  }
+};
+
+#endif // VORONOI_SOLVER_HPP