True Adding second order in space for 3 current solver with induced limitation

src/scheme/CellbyCellLimitation.hpp

+#ifndef CELL_BY_CELL_LIMITATION_HPP
+#define CELL_BY_CELL_LIMITATION_HPP
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+
+#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshEdgeBoundary.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatEdgeBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/PugsTraits.hpp>
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class CellByCellLimitation
+{
+ private:
+  using MeshType                    = MeshTypeT;
+  static constexpr size_t Dimension = MeshType::Dimension;
+
+  using Rdxd = TinyMatrix<Dimension>;
+  using Rd   = TinyVector<Dimension>;
+
+ public:
+  void
+  computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
+                                                            const CellValue<double>& q,
+                                                            const CellValue<Rd>& gradq,
+                                                            CellValue<double>& Limitor_q,
+                                                            const bool enableWeakBoundPositivityOnly = false) const
+  {
+    Assert(Dimension == 2 or Dimension == 3);
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+
+    const auto& xr = mesh.xr();
+
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    // auto volumes_cell = mesh_data.Vj();
+
+    const auto& xcentroid = mesh_data.xj();
+    const auto& xface     = mesh_data.xl();
+
+    if (Dimension == 2) {
+      parallel_for(
+        mesh.numberOfCells(),
+        PUGS_LAMBDA(CellId j) {   // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+          const double State = q[j];
+          const Rd& Cent     = xcentroid[j];   // mesh.xi(j);
+          const Rd& Gradl    = gradq[j];
+
+          // Min et Max des valeurs voisines
+
+          double minvalvois(State);
+          double maxvalvois(State);
+
+          // Calcul des extremas de la fonction
+
+          double minval(std::numeric_limits<double>::max());
+          double maxval(-std::numeric_limits<double>::max());
+
+          // At nodes
+          const auto& cell_to_node = cell_to_node_matrix[j];
+          // At faces
+          const auto& cell_to_face = cell_to_face_matrix[j];
+
+          for (size_t l = 0; l < cell_to_node.size(); ++l) {
+            const NodeId& node       = cell_to_node[l];
+            const auto& node_to_cell = node_to_cell_matrix[node];
+
+            for (size_t k = 0; k < node_to_cell.size(); ++k) {
+              const CellId& cellvois     = node_to_cell[k];
+              const double statevoiscell = q[cellvois];
+              const Rd Centvois          = xcentroid[cellvois];   // mesh.element()[mayvois].centroid().x();
+              const Rd Mvois(Centvois - Cent);
+
+              const double valVois = State + dot(Gradl, Mvois);   // Partie MinMod
+
+              minval = std::min(minval, valVois);
+              maxval = std::max(maxval, valVois);
+
+              minvalvois = std::min(minvalvois, statevoiscell);
+              maxvalvois = std::max(maxvalvois, statevoiscell);
+            }
+
+            const Rd M(xr[node] - Cent);
+
+            const double valLineO2node = State + dot(Gradl, M);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2node);
+            maxval = std::max(maxval, valLineO2node);
+          }
+
+          for (size_t l = 0; l < cell_to_face.size(); ++l) {
+            const FaceId& face = cell_to_face[l];
+
+            const Rd Mb(xface[face] - Cent);
+
+            const double valLineO2bras = State + dot(Gradl, Mb);   // Idem aux demi faces
+
+            minval = std::min(minval, valLineO2bras);
+            maxval = std::max(maxval, valLineO2bras);
+          }
+
+          static const double epsZer0 = 1e-15;
+          // on applique le traitement idem de l'ordre 2..
+          double coef1 = 0, coef2 = 0;
+          if (enableWeakBoundPositivityOnly) {
+            coef1                  = 1;
+            const double minGlobal = 1e-12;
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minGlobal - State) / ((minval - State));
+
+          } else {
+            if (std::fabs(maxval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef1 = 1.;
+            else
+              coef1 = (maxvalvois - State) / ((maxval - State));
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minvalvois - State) / ((minval - State));
+          }
+          const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+          Limitor_q[j] = alfa;
+          // grad(i) *= alfa;
+        });
+    } else {
+      // throw NormalError(" Limiteur Scal Grad non encore implemente 3D");
+      const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+      const auto& xedge               = mesh_data.xe();
+
+      parallel_for(
+        mesh.numberOfCells(),
+        PUGS_LAMBDA(CellId j) {   // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+          const double State = q[j];
+          const Rd& Cent     = xcentroid[j];   // mesh.xi(j);
+          const Rd& Gradl    = gradq[j];
+
+          // Min et Max des valeurs voisines
+
+          double minvalvois(State);
+          double maxvalvois(State);
+
+          // Calcul des extremas de la fonction
+
+          double minval(std::numeric_limits<double>::max());
+          double maxval(-std::numeric_limits<double>::max());
+
+          // At nodes
+          const auto& cell_to_node = cell_to_node_matrix[j];
+          // At faces
+          const auto& cell_to_face = cell_to_face_matrix[j];
+          // At faces
+          const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+          for (size_t l = 0; l < cell_to_node.size(); ++l) {
+            const NodeId& node       = cell_to_node[l];
+            const auto& node_to_cell = node_to_cell_matrix[node];
+
+            for (size_t k = 0; k < node_to_cell.size(); ++k) {
+              const CellId& cellvois     = node_to_cell[k];
+              const double statevoiscell = q[cellvois];
+              const Rd Centvois          = xcentroid[cellvois];   // mesh.element()[mayvois].centroid().x();
+              const Rd Mvois(Centvois - Cent);
+
+              const double valVois = State + dot(Gradl, Mvois);   // Partie MinMod
+
+              minval = std::min(minval, valVois);
+              maxval = std::max(maxval, valVois);
+
+              minvalvois = std::min(minvalvois, statevoiscell);
+              maxvalvois = std::max(maxvalvois, statevoiscell);
+            }
+
+            const Rd M(xr[node] - Cent);
+
+            const double valLineO2node = State + dot(Gradl, M);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2node);
+            maxval = std::max(maxval, valLineO2node);
+          }
+
+          for (size_t l = 0; l < cell_to_face.size(); ++l) {
+            const FaceId& face = cell_to_face[l];
+
+            const Rd Mb(xface[face] - Cent);
+
+            const double valLineO2face = State + dot(Gradl, Mb);   // Idem aux demi faces
+
+            minval = std::min(minval, valLineO2face);
+            maxval = std::max(maxval, valLineO2face);
+          }
+
+          for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+            const EdgeId& edge = cell_to_edge[l];
+
+            const Rd Me(xedge[edge] - Cent);
+
+            const double valLineO2edge = State + dot(Gradl, Me);   // Idem aux demi faces
+
+            minval = std::min(minval, valLineO2edge);
+            maxval = std::max(maxval, valLineO2edge);
+          }
+
+          static const double epsZer0 = 1e-15;
+          // on applique le traitement idem de l'ordre 2..
+          double coef1 = 0, coef2 = 0;
+          if (enableWeakBoundPositivityOnly) {
+            coef1                  = 1;
+            const double minGlobal = 1e-12;
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minGlobal - State) / ((minval - State));
+
+          } else {
+            if (std::fabs(maxval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef1 = 1.;
+            else
+              coef1 = (maxvalvois - State) / ((maxval - State));
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minvalvois - State) / ((minval - State));
+          }
+          const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+          Limitor_q[j] = alfa;
+          // grad(i) *= alfa;
+        });
+    }
+  }
+
+  void
+  computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
+                                                                const CellValue<double>& rho,
+                                                                const CellValue<Rd>& gradrho,
+                                                                const CellValue<double>& internalnrj,
+                                                                const CellValue<Rd>& gradinternalnrj,
+                                                                // const CellValue<Rd>& u,
+                                                                const CellValue<Rdxd>& gradu,
+                                                                CellValue<double>& Limitor_eps,
+                                                                const bool enableWeakBoundPositivityOnly = false) const
+  {
+    Assert(Dimension == 2 or Dimension == 3);
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+    const auto& xr                  = mesh.xr();
+
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const auto& xcentroid = mesh_data.xj();
+    const auto& xface     = mesh_data.xl();
+
+    if (Dimension == 2) {
+      parallel_for(
+        mesh.numberOfCells(),
+        PUGS_LAMBDA(CellId j) {   // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+          const double rhol   = rho[j];
+          const Rd& Cent      = xcentroid[j];
+          const Rd& grad_rhol = gradrho[j];
+
+          // const Rd& ul        = u[j];
+          const Rdxd& grad_ul = gradu[j];
+
+          const double internal_nrj = internalnrj[j];
+          // const double kinetic_nrj  = .5 * dot(ul, ul);
+
+          const Rd& grad_internal_nrj = gradinternalnrj[j];
+
+          // const Rd& grad_kinetic_internal_nrj = gradeps[j];
+
+          // const double qtt_limitor_density = rhol / (rhol + dot(grad,(x-xc));
+          // At nodes
+          const auto& cell_to_node = cell_to_node_matrix[j];
+          // At faces
+          const auto& cell_to_face = cell_to_face_matrix[j];
+
+          // Min et Max des valeurs voisines
+          const double State = internal_nrj;
+
+          double minvalvois(State);
+          double maxvalvois(State);
+
+          // Calcul des extremas de la fonction
+
+          double minval(std::numeric_limits<double>::max());
+          double maxval(-std::numeric_limits<double>::max());
+
+          for (size_t l = 0; l < cell_to_node.size(); ++l) {
+            const NodeId& node       = cell_to_node[l];
+            const auto& node_to_cell = node_to_cell_matrix[node];
+
+            for (size_t k = 0; k < node_to_cell.size(); ++k) {
+              const CellId& cellvois = node_to_cell[k];
+              const Rd Centvois      = xcentroid[cellvois];
+              const Rd Mvois(Centvois - Cent);
+
+              const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mvois));
+
+              const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Mvois);
+              const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mvois);
+
+              const double valVois = internal_nrj + High_orderinternalnrj -
+                                     .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie MinMod
+
+              minval = std::min(minval, valVois);
+              maxval = std::max(maxval, valVois);
+
+              const double statevoiscell = internalnrj[cellvois];
+
+              minvalvois = std::min(minvalvois, statevoiscell);
+              maxvalvois = std::max(maxvalvois, statevoiscell);
+            }
+
+            const Rd Mr(xr[node] - Cent);
+
+            const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mr));
+
+            const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Mr);
+            const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mr);
+
+            const double valLineO2node = internal_nrj + High_orderinternalnrj -
+                                         .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2node);
+            maxval = std::max(maxval, valLineO2node);
+          }
+
+          for (size_t l = 0; l < cell_to_face.size(); ++l) {
+            const FaceId& face = cell_to_face[l];
+
+            const Rd Mb(xface[face] - Cent);
+
+            const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mb));
+
+            const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Mb);
+            const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mb);
+
+            const double valLineO2face = internal_nrj + High_orderinternalnrj -
+                                         .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2face);
+            maxval = std::max(maxval, valLineO2face);
+          }
+
+          static const double epsZer0 = 1e-15;
+
+          // on applique le traitement idem de l'ordre 2..
+          double coef1 = 0, coef2 = 0;
+          if (enableWeakBoundPositivityOnly) {
+            coef1                  = 1;
+            const double minGlobal = 1e-12;
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minGlobal - State) / ((minval - State));
+
+          } else {
+            if (std::fabs(maxval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef1 = 1.;
+            else
+              coef1 = (maxvalvois - State) / ((maxval - State));
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minvalvois - State) / ((minval - State));
+          }
+          const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+          Limitor_eps[j] = alfa;
+        });
+    } else {
+      // throw NormalError(" Limiteur Scal Grad internal nrj non encore implemente 3D");
+      const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+
+      const auto& xedge = mesh_data.xe();
+
+      parallel_for(
+        mesh.numberOfCells(),
+        PUGS_LAMBDA(CellId j) {   // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+          const double rhol   = rho[j];
+          const Rd& Cent      = xcentroid[j];
+          const Rd& grad_rhol = gradrho[j];
+
+          // const Rd& ul        = u[j];
+          const Rdxd& grad_ul = gradu[j];
+
+          const double internal_nrj = internalnrj[j];
+          // const double kinetic_nrj  = .5 * dot(ul, ul);
+
+          const Rd& grad_internal_nrj = gradinternalnrj[j];
+
+          // const Rd& grad_kinetic_internal_nrj = gradeps[j];
+
+          // const double qtt_limitor_density = rhol / (rhol + dot(grad,(x-xc));
+          // At nodes
+          const auto& cell_to_node = cell_to_node_matrix[j];
+          // At faces
+          const auto& cell_to_face = cell_to_face_matrix[j];
+          // At edges
+          const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+          // Min et Max des valeurs voisines
+          const double State = internal_nrj;
+
+          double minvalvois(State);
+          double maxvalvois(State);
+
+          // Calcul des extremas de la fonction
+
+          double minval(std::numeric_limits<double>::max());
+          double maxval(-std::numeric_limits<double>::max());
+
+          for (size_t l = 0; l < cell_to_node.size(); ++l) {
+            const NodeId& node       = cell_to_node[l];
+            const auto& node_to_cell = node_to_cell_matrix[node];
+
+            for (size_t k = 0; k < node_to_cell.size(); ++k) {
+              const CellId& cellvois = node_to_cell[k];
+              const Rd Centvois      = xcentroid[cellvois];
+              const Rd Mvois(Centvois - Cent);
+
+              const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mvois));
+
+              const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Mvois);
+              const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mvois);
+
+              const double valVois = internal_nrj + High_orderinternalnrj -
+                                     .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie MinMod
+
+              minval = std::min(minval, valVois);
+              maxval = std::max(maxval, valVois);
+
+              const double statevoiscell = internalnrj[cellvois];
+
+              minvalvois = std::min(minvalvois, statevoiscell);
+              maxvalvois = std::max(maxvalvois, statevoiscell);
+            }
+
+            const Rd Mr(xr[node] - Cent);
+
+            const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mr));
+
+            const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Mr);
+            const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mr);
+
+            const double valLineO2node = internal_nrj + High_orderinternalnrj -
+                                         .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2node);
+            maxval = std::max(maxval, valLineO2node);
+          }
+
+          for (size_t l = 0; l < cell_to_face.size(); ++l) {
+            const FaceId& face = cell_to_face[l];
+
+            const Rd Mb(xface[face] - Cent);
+
+            const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mb));
+
+            const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Mb);
+            const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mb);
+
+            const double valLineO2face = internal_nrj + High_orderinternalnrj -
+                                         .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2face);
+            maxval = std::max(maxval, valLineO2face);
+          }
+
+          for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+            const EdgeId& edge = cell_to_edge[l];
+
+            const Rd Me(xedge[edge] - Cent);
+
+            const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Me));
+
+            const Rd& High_ordervelocity       = qtt_limitor_density * (grad_ul * Me);
+            const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Me);
+
+            const double valLineO2edge = internal_nrj + High_orderinternalnrj -
+                                         .5 * dot(High_ordervelocity, High_ordervelocity);   // Partie Dukowicz
+
+            minval = std::min(minval, valLineO2edge);
+            maxval = std::max(maxval, valLineO2edge);
+          }
+
+          static const double epsZer0 = 1e-15;
+
+          // on applique le traitement idem de l'ordre 2..
+          double coef1 = 0, coef2 = 0;
+          if (enableWeakBoundPositivityOnly) {
+            coef1                  = 1;
+            const double minGlobal = 1e-12;
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minGlobal - State) / ((minval - State));
+
+          } else {
+            if (std::fabs(maxval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef1 = 1.;
+            else
+              coef1 = (maxvalvois - State) / ((maxval - State));
+
+            if (std::fabs(minval - State) <= epsZer0)   // epsIsEqualAbs())
+              coef2 = 1.;
+            else
+              coef2 = (minvalvois - State) / ((minval - State));
+          }
+          const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+          Limitor_eps[j] = alfa;
+        });
+    }
+  }
+};
+#endif

src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.hpp

+#ifndef ROE_VISCOUS_FORM_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+#define ROE_VISCOUS_FORM_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/CellbyCellLimitation.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+#include <memory>
+#include <vector>
+
+// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+//                  const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+roeViscousFormEulerianCompositeSolver_v2_o2(
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+  const double& gamma,
+  const std::shared_ptr<const DiscreteFunctionVariant>& c,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+  const size_t& degree,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const double& dt,
+  const bool check = false);
+
+#endif   // ROE_VISCOUS_FORM_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP

src/scheme/RusanovEulerianCompositeSolver_o2.hpp

+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_O2_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_O2_HPP
+
+#include <memory>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/CellbyCellLimitation.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+#include <vector>
+
+// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+//                  const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver_o2(
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+  const double& gamma,
+  const std::shared_ptr<const DiscreteFunctionVariant>& c,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+  const size_t& degree,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const double& dt,
+  const bool check = false);
+
+#endif   // RUSANOV_EULERIAN_COMPOSITE_SOLVER_HPP

src/scheme/RusanovEulerianCompositeSolver_v2_o2.hpp

+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+
+#include <memory>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/CellbyCellLimitation.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+#include <vector>
+
+// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+//                  const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver_v2_o2(
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+  const std::shared_ptr<const DiscreteFunctionVariant>& u,
+  const std::shared_ptr<const DiscreteFunctionVariant>& E,
+  const double& gamma,
+  const std::shared_ptr<const DiscreteFunctionVariant>& c,
+  const std::shared_ptr<const DiscreteFunctionVariant>& p,
+  const size_t& degree,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const double& dt,
+  const bool check = false);
+
+#endif   // RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_O2_xHPP