Add entropic_eucclhyd_solver function

The code works in 1D/2D but it is very expensive. Probably due to the
allocations in the Polynomial class

CMakeLists.txt

@@ -545,7 +545,6 @@ target_link_libraries(
   pugs
   PugsMesh
   PugsAlgebra
-  PugsAnalysis
   PugsUtils
   PugsLanguage
   PugsLanguageAST
@@ -555,6 +554,7 @@ target_link_libraries(
   PugsAlgebra
   PugsAnalysis
   PugsScheme
+  PugsAnalysis
   PugsUtils
   PugsOutput
   PugsLanguageUtils

src/analysis/Polynomial1D.hpp

@@ -263,6 +263,15 @@ class [[nodiscard]] Polynomial1D
     for (size_t i = 0; i < coeffients.size(); ++i) {
       m_coefficients[i] = coeffients[i];
     }
+
+    size_t real_degree = this->_getRealDegree();
+    if (real_degree + 1 < m_coefficients.size()) {
+      Array<double> real_coefficients(real_degree + 1);
+      for (size_t i = 0; i <= real_degree; ++i) {
+        real_coefficients[i] = m_coefficients[i];
+      }
+      m_coefficients = real_coefficients;
+    }
   }
 
   PUGS_INLINE

src/language/modules/SchemeModule.cpp

@@ -188,6 +188,42 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("entropic_eucclhyd_solver",
+                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
+                              double, std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::vector<std::shared_ptr<
+                                                                          const IBoundaryConditionDescriptor>>&,
+                                                                        const double&)>>(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                 const std::shared_ptr<const IDiscreteFunction>& u,
+                                 const std::shared_ptr<const IDiscreteFunction>& E,
+                                 const std::shared_ptr<const IDiscreteFunction>& c,
+                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt_max)
+                                -> std::tuple<
+                                  double, std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+                                  std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> {
+                                return AcousticSolverHandler{AcousticSolverHandler::SolverType::Eucclhyd,
+                                                             rho,
+                                                             c,
+                                                             u,
+                                                             p,
+                                                             bc_descriptor_list}
+                                  .solver()
+                                  .apply_computing_entropic_dt(dt_max, rho, u, E, p);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("eucclhyd_solver",
                             std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                               std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,

src/scheme/AcousticSolver.cpp

 #include <scheme/AcousticSolver.hpp>
 
+#include <analysis/Polynomial1D.hpp>
+#include <analysis/Polynomial1DRootComputer.hpp>
 #include <language/utils/InterpolateItemValue.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
@@ -84,6 +86,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   NodeValue<const Rd> m_ur;
   NodeValuePerCell<const Rd> m_Fjr;
 
+  NodeValuePerCell<const Rdxd> m_Ajr;
+
   CellValue<const double>
   _getRhoC(const DiscreteScalarFunction& rho, const DiscreteScalarFunction& c)
   {
@@ -393,6 +397,294 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     return F;
   }
 
+  double
+  _computeMaxDensityDt(const MeshType& mesh,
+                       const DiscreteFunctionP0<Dimension, double>& rho,
+                       const double dt_max) const
+  {
+    double dt = dt_max;
+
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const CellValue<const double> Vj     = mesh_data.Vj();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    CellValue<double> dVj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double dV = 0;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          dV += dot(Cjr(j, R), m_ur[r]);
+        }
+
+        dVj[j] = dV;
+      });
+
+    CellValue<double> Gj{mesh.connectivity()};
+
+    if constexpr (Dimension == 2) {
+      for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double G = 0;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r   = cell_nodes[R];
+          const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+          G -= dot(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]});
+        }
+        Gj[j] = 0.5 * G;
+      }
+    }
+
+    for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+      Polynomial1D p_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])});
+      if constexpr (Dimension == 2) {
+        p_tau = p_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])});
+      }
+      Polynomial1DRootComputer computer{p_tau, 0, dt};
+      std::optional dt_tau = computer.getFirstRoot();
+      if (dt_tau.has_value()) {
+        Assert(dt_tau.value() <= dt);
+        dt = dt_tau.value();
+      }
+    }
+
+    if (dt != dt_max) {
+      std::cout << "volume variation imposes time step\n";
+    }
+
+    if (dt < dt_max) {
+      dt *= 0.95 * dt;
+    }
+
+    return parallel::allReduceMin(dt);
+  }
+
+  double
+  _computeMaxEnergyDt(const MeshType& mesh,
+                      const DiscreteScalarFunction& rho,
+                      const DiscreteVectorFunction& u,
+                      const DiscreteScalarFunction& E,
+                      const DiscreteScalarFunction& p,
+                      const double dt_max) const
+  {
+    double dt = dt_max;
+
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const CellValue<const double> Vj     = mesh_data.Vj();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    CellValue<double> dVj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double dV = 0;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          dV += dot(Cjr(j, R), m_ur[r]);
+        }
+
+        dVj[j] = dV;
+      });
+
+    CellValue<double> dSj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double dS = 0;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+
+          Rd du = u[j] - m_ur[r];
+          dS += dot(m_Ajr(j, R) * du, du);
+        }
+
+        dSj[j] = dS;
+      });
+
+    CellValue<Rd> dPj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd dP = zero;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+
+          Rd du = u[j] - m_ur[r];
+          dP += m_Ajr(j, R) * du;
+        }
+
+        dPj[j] = dP;
+      });
+
+    for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+      const double inv_mj = 1 / (rho[j] * Vj[j]);
+
+      Polynomial1D p_e(
+        {E[j] - 0.5 * dot(u[j], u[j]), inv_mj * (dSj[j] - p[j] * dVj[j]), -inv_mj * inv_mj * dot(dPj[j], dPj[j])});
+      Polynomial1DRootComputer computer{p_e, 0, dt};
+      std::optional dt_e = computer.getFirstRoot();
+      if (dt_e.has_value()) {
+        Assert(dt_e.value() <= dt);
+        dt = dt_e.value();
+      }
+    }
+
+    if (dt < dt_max) {
+      dt *= 0.95 * dt;
+    }
+
+    if (dt != dt_max) {
+      std::cout << "internal energy variation imposes time step\n";
+    }
+
+    return parallel::allReduceMin(dt);
+  }
+
+  double
+  _computeMaxEntropyDt(const MeshType& mesh,
+                       const DiscreteScalarFunction& rho,
+                       const DiscreteVectorFunction& u,
+                       const DiscreteScalarFunction& p,
+                       const double dt_max) const
+  {
+    double dt = dt_max;
+
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const CellValue<const double> Vj     = mesh_data.Vj();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    CellValue<double> dVj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double dV = 0;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+          dV += dot(Cjr(j, R), m_ur[r]);
+        }
+
+        dVj[j] = dV;
+      });
+
+    CellValue<double> Gj{mesh.connectivity()};
+
+    if constexpr (Dimension == 2) {
+      for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double G = 0;
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r   = cell_nodes[R];
+          const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+          G -= dot(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]});
+        }
+        Gj[j] = 0.5 * G;
+      }
+    }
+
+    CellValue<double> dSj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        double dS = 0;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+
+          const Rd du = u[j] - m_ur[r];
+          dS += dot(m_Ajr(j, R) * du, du);
+        }
+
+        dSj[j] = dS;
+      });
+
+    CellValue<Rd> dPj{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        Rd dP = zero;
+
+        for (size_t R = 0; R < cell_nodes.size(); ++R) {
+          const NodeId r = cell_nodes[R];
+
+          Rd du = u[j] - m_ur[r];
+          dP += m_Ajr(j, R) * du;
+        }
+
+        dPj[j] = dP;
+      });
+
+#warning fixed gamma value
+    const double gamma = 1.4;
+
+    for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+      const double inv_mj = 1 / (rho[j] * Vj[j]);
+      const double inv_Vj = 1 / Vj[j];
+
+      Polynomial1D DT({0, 1});
+      std::vector<double> delta_S = {dSj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])};
+      if constexpr (Dimension == 2) {
+        delta_S[1] += p[j] * Gj[j];
+      }
+      Polynomial1D p_S0(delta_S);
+
+      std::vector<double> delta_V{dVj[j] /*, dG */};
+      if constexpr (Dimension == 2) {
+        delta_V.push_back(Gj[j]);
+      }
+      Polynomial1D p_deltaV(delta_V);
+
+      Polynomial1D p_S1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])});
+
+      Polynomial1D p_S =                                                                                 //
+        p_S0                                                                                             //
+        + (inv_Vj * DT) * (((gamma - 1) * p_S1) * p_deltaV + (gamma - 2) * p[j] * p_deltaV * p_deltaV)   //
+        + (inv_Vj * DT) * (inv_Vj * DT) * (((gamma - 1) * (gamma - 2)) * p_S1) * p_deltaV * p_deltaV;
+
+      Polynomial1DRootComputer computer{p_S, 0, dt};
+      std::optional dt_S = computer.getFirstRoot();
+      if (dt_S.has_value()) {
+        Assert(dt_S.value() <= dt);
+        dt = dt_S.value();
+      }
+    }
+
+    if (dt != dt_max) {
+      std::cout << "entropy variation imposes time step\n";
+    }
+    return parallel::allReduceMin(dt);
+  }
+
   AcousticSolver(SolverType solver_type,
                  const std::shared_ptr<const MeshType>& p_mesh,
                  const DiscreteScalarFunction& rho,
@@ -407,6 +699,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     CellValue<const double> rhoc     = this->_getRhoC(rho, c);
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rhoc);
 
+    m_Ajr = Ajr;
+
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
     NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
 
@@ -494,6 +788,91 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                        *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
   }
 
+  std::tuple<double,
+             std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteScalarFunction>,
+             std::shared_ptr<const DiscreteVectorFunction>,
+             std::shared_ptr<const DiscreteScalarFunction>>
+  apply_computing_entropic_dt(const double& dt_max,
+                              const std::shared_ptr<const MeshType>& mesh,
+                              const DiscreteScalarFunction& rho,
+                              const DiscreteVectorFunction& u,
+                              const DiscreteScalarFunction& E,
+                              const DiscreteScalarFunction& p) const
+  {
+    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
+    double dt = dt_max;
+
+    dt = _computeMaxDensityDt(*mesh, rho, dt);
+    dt = _computeMaxEnergyDt(*mesh, rho, u, E, p, dt);
+    dt = _computeMaxEntropyDt(*mesh, rho, u, p, dt);
+
+    NodeValue<Rd> new_xr = copy(mesh->xr());
+    parallel_for(
+      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(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 += m_Fjr(j, R);
+          energy_fluxes += dot(m_Fjr(j, R), m_ur[r]);
+        }
+        const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+        new_u[j] -= dt_over_Mj * momentum_fluxes;
+        new_E[j] -= dt_over_Mj * energy_fluxes;
+      });
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {dt, new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
+            std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
+            std::make_shared<DiscreteScalarFunction>(new_mesh, new_E)};
+  }
+
+  std::tuple<double,
+             std::shared_ptr<const IMesh>,
+             std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>>
+  apply_computing_entropic_dt(const double& dt_max,
+                              const std::shared_ptr<const IDiscreteFunction>& rho,
+                              const std::shared_ptr<const IDiscreteFunction>& u,
+                              const std::shared_ptr<const IDiscreteFunction>& E,
+                              const std::shared_ptr<const IDiscreteFunction>& p) const
+  {
+    std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
+    if (not i_mesh) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    return this->apply_computing_entropic_dt(dt_max, std::dynamic_pointer_cast<const MeshType>(i_mesh),
+                                             *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
+                                             *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
+                                             *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E),
+                                             *std::dynamic_pointer_cast<const DiscreteScalarFunction>(p));
+  }
+
   AcousticSolver(SolverType solver_type,
                  const std::shared_ptr<const IMesh>& mesh,
                  const std::shared_ptr<const IDiscreteFunction>& rho,

src/scheme/AcousticSolver.hpp

@@ -32,6 +32,17 @@ class AcousticSolverHandler
           const std::shared_ptr<const IDiscreteFunction>& u,
           const std::shared_ptr<const IDiscreteFunction>& E) const = 0;
 
+    virtual std::tuple<double,
+                       std::shared_ptr<const IMesh>,
+                       std::shared_ptr<const IDiscreteFunction>,
+                       std::shared_ptr<const IDiscreteFunction>,
+                       std::shared_ptr<const IDiscreteFunction>>
+    apply_computing_entropic_dt(const double& dt_max,
+                                const std::shared_ptr<const IDiscreteFunction>& rho,
+                                const std::shared_ptr<const IDiscreteFunction>& u,
+                                const std::shared_ptr<const IDiscreteFunction>& E,
+                                const std::shared_ptr<const IDiscreteFunction>& p) const = 0;
+
     IAcousticSolver()                  = default;
     IAcousticSolver(IAcousticSolver&&) = default;
     IAcousticSolver& operator=(IAcousticSolver&&) = default;

src/scheme/CMakeLists.txt

@@ -6,3 +6,9 @@ add_library(
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
   DiscreteFunctionVectorInterpoler.cpp)
+
+# Additional dependencies
+add_dependencies(
+  PugsScheme
+  PugsAnalysis
+)

tests/CMakeLists.txt

@@ -130,8 +130,8 @@ target_link_libraries (unit_tests
   PugsLanguage
   PugsMesh
   PugsAlgebra
-  PugsAnalysis
   PugsScheme
+  PugsAnalysis
   PugsOutput
   PugsUtils
   kokkos
@@ -153,10 +153,10 @@ target_link_libraries (mpi_unit_tests
   PugsLanguageAlgorithms
   PugsMesh
   PugsAlgebra
-  PugsAnalysis
   PugsUtils
   PugsLanguageUtils
   PugsScheme
+  PugsAnalysis
   PugsOutput
   PugsUtils
   PugsAlgebra