[ci-skip] New files for 2D kinetic solver

src/language/modules/SchemeModule.cpp

@@ -33,6 +33,7 @@
 #include <scheme/DiscreteFunctionVectorIntegrator.hpp>
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
 #include <scheme/EulerKineticSolver.hpp>
+#include <scheme/EulerKineticSolver2D.hpp>
 #include <scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.hpp>
 #include <scheme/EulerKineticSolverAcousticLagrangeFV.hpp>
 #include <scheme/EulerKineticSolverFirstOrderFV.hpp>
@@ -930,6 +931,39 @@ SchemeModule::SchemeModule()
 
                                           ));
 
+  this->_addBuiltinFunction("euler_kinetic_2D",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector,
+                                 const double& eps, const size_t& SpaceOrder,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolver2D(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho, F_rho_u1,
+                                                            F_rho_u2, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_2D",
+                            std::function(
+
+                              [](const std::vector<double>& lambda_vector,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWaves2D(lambda_vector, rho, rho_u, rho_E, gamma);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("euler_kinetic_first_order_FV",
                             std::function(
 

src/scheme/CMakeLists.txt

@@ -22,6 +22,7 @@ add_library(
   EulerKineticSolverOneFluxMood.cpp
   EulerKineticSolverMoodFD.cpp
   EulerKineticSolverMoodFV.cpp
+  EulerKineticSolver2D.cpp
   EulerKineticSolverThirdOrderMoodFV.cpp
   EulerKineticSolverThirdOrderMoodFD.cpp
   EulerKineticSolverThirdOrderFV.cpp

src/scheme/EulerKineticSolver2D.cpp

+#include <scheme/EulerKineticSolver2D.hpp>
+
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <tuple>
+
+template <MeshConcept MeshType>
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+getEulerKineticWaves2D(const std::shared_ptr<const MeshType> p_mesh,
+                       const std::vector<double> lambda_vector,
+                       DiscreteFunctionP0<const double> rho,
+                       DiscreteFunctionP0<const double> rho_u,
+                       DiscreteFunctionP0<const double> rho_E,
+                       const double gamma)
+{
+  if (lambda_vector.size() < 2) {
+    throw NormalError("lambda vector must be of dimension greater than 2");
+  }
+
+  MeshData<MeshType>& mesh_data                               = MeshDataManager::instance().getMeshData(*p_mesh);
+  const NodeValue<const TinyVector<MeshType::Dimension>> m_xr = p_mesh->xr();
+  const CellValue<const TinyVector<MeshType::Dimension>> m_xj = mesh_data.xj();
+  const CellValue<const double> m_dx_table                    = mesh_data.Vj();
+
+  DiscreteFunctionP0Vector<double> Fn_rho(p_mesh, lambda_vector.size());
+  DiscreteFunctionP0Vector<double> Fn_rho_u(p_mesh, lambda_vector.size());
+  DiscreteFunctionP0Vector<double> Fn_rho_E(p_mesh, lambda_vector.size());
+
+  const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+  CellArray<double> vec_A{p_mesh->connectivity(), 3};
+  CellValue<double> p{p_mesh->connectivity()};
+
+  for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+    p[cell_id]        = (gamma - 1) * (rho_E[cell_id] - 0.5 * rho_u[cell_id] * rho_u[cell_id] / rho[cell_id]);
+    vec_A[cell_id][0] = rho_u[cell_id];
+    vec_A[cell_id][1] = rho_u[cell_id] * rho_u[cell_id] / rho[cell_id] + p[cell_id];
+    vec_A[cell_id][2] = (rho_E[cell_id] + p[cell_id]) * rho_u[cell_id] / rho[cell_id];
+  }
+
+  const size_t nb_waves           = lambda_vector.size();
+  double inv_L2norm_lambda_vector = 0;
+  for (size_t i = 0; i < nb_waves; ++i) {
+    inv_L2norm_lambda_vector += std::pow(lambda_vector[i], 2);
+  }
+  inv_L2norm_lambda_vector = 1. / inv_L2norm_lambda_vector;
+
+  for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+    for (size_t i = 0; i < nb_waves; ++i) {
+      Fn_rho[cell_id][i] =
+        (1. / nb_waves) * rho[cell_id] + inv_L2norm_lambda_vector * vec_A[cell_id][0] * lambda_vector[i];
+      Fn_rho_u[cell_id][i] =
+        (1. / nb_waves) * rho_u[cell_id] + inv_L2norm_lambda_vector * vec_A[cell_id][1] * lambda_vector[i];
+      Fn_rho_E[cell_id][i] =
+        (1. / nb_waves) * rho_E[cell_id] + inv_L2norm_lambda_vector * vec_A[cell_id][2] * lambda_vector[i];
+    }
+  }
+
+  return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(Fn_rho),
+                         std::make_shared<DiscreteFunctionVariant>(Fn_rho_u),
+                         std::make_shared<DiscreteFunctionVariant>(Fn_rho_E));
+}
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+getEulerKineticWaves2D(const std::vector<double>& lambda_vector,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho_u,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
+                       const double& gamma)
+{
+  const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({rho, rho_u, rho_E});
+  if (mesh_v.use_count() == 0) {
+    throw NormalError("rho rho_u and rho_E have to be defined on the same mesh");
+  }
+
+  return std::visit(
+    [&](auto&& p_mesh)
+      -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>> {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 1)) {
+        auto rho_h   = rho->get<DiscreteFunctionP0<const double>>();
+        auto rho_u_h = rho_u->get<DiscreteFunctionP0<const double>>();
+        auto rho_E_h = rho_E->get<DiscreteFunctionP0<const double>>();
+
+        return getEulerKineticWaves2D(p_mesh, lambda_vector, rho_h, rho_u_h, rho_E_h, gamma);
+      } else {
+        throw NotImplementedError("Invalid mesh type for Scalar Kinetic solver");
+      }
+    },
+    mesh_v->variant());
+}
+
+template <MeshConcept MeshType>
+class EulerKineticSolver2D
+{
+ private:
+  constexpr static size_t Dimension = MeshType::Dimension;
+  using Rd                          = TinyVector<Dimension>;
+
+  const double m_dt;
+  const double m_eps;
+  const std::vector<TinyVector<Dimension>>& m_lambda_vector;
+  const size_t m_SpaceOrder;
+  const DiscreteFunctionP0Vector<const double> m_Fn_rho;
+  const DiscreteFunctionP0Vector<const double> m_Fn_rho_u1;
+  const DiscreteFunctionP0Vector<const double> m_Fn_rho_u2;
+  const DiscreteFunctionP0Vector<const double> m_Fn_rho_E;
+  const double m_gamma;
+  const std::shared_ptr<const MeshType> p_mesh;
+  const MeshType& m_mesh;
+  const CellValue<const double> m_dx_table = MeshDataManager::instance().getMeshData(m_mesh).Vj();
+  CellValue<const double> m_inv_dx_table;
+  const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+  const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+
+ public:
+  CellArray<TinyVector<MeshType::Dimension>>
+  getA(const CellValue<const double>& rho,
+       const CellValue<const TinyVector<Dimension>>& rho_u,
+       const CellValue<const double>& rho_E) const
+  {
+    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+    CellArray<TinyVector<Dimension>> vec_A{m_mesh.connectivity(), 4};
+    CellValue<double> p{m_mesh.connectivity()};
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        vec_A[cell_id][0][0] = rho_u[cell_id][0];
+        vec_A[cell_id][1][0] = rho_u[cell_id][0] * rho_u[cell_id][0] / rho[cell_id] + p[cell_id];
+        vec_A[cell_id][2][0] = rho_u[cell_id][0] * rho_u[cell_id][1] / rho[cell_id];
+        vec_A[cell_id][3][0] = (rho_E[cell_id] + p[cell_id]) * rho_u[cell_id][0] / rho[cell_id];
+
+        vec_A[cell_id][0][1] = rho_u[cell_id][1];
+        vec_A[cell_id][1][1] = rho_u[cell_id][0] * rho_u[cell_id][1] / rho[cell_id];
+        vec_A[cell_id][2][1] = rho_u[cell_id][1] * rho_u[cell_id][1] / rho[cell_id] + p[cell_id];
+        vec_A[cell_id][3][1] = (rho_E[cell_id] + p[cell_id]) * rho_u[cell_id][1] / rho[cell_id];
+      });
+    return vec_A;
+  }
+
+  CellArray<double>
+  compute_M(const CellValue<const double>& u, const CellValue<TinyVector<MeshType::Dimension>> A_u)
+  {
+    const size_t nb_waves = m_lambda_vector.size();
+    CellArray<double> M{m_mesh.connectivity(), nb_waves};
+
+    TinyVector<MeshType::Dimension> inv_S = zero;
+    for (size_t i = 0; i < nb_waves; ++i) {
+      for (size_t d = 0; d < MeshType::Dimension; ++d) {
+        inv_S[d] += std::pow(m_lambda_vector[i][d], 2);
+      }
+    }
+    inv_S = 1. / inv_S;
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        for (size_t i = 0; i < nb_waves; ++i) {
+          M[cell_id][i] = (1. / nb_waves) * u[cell_id] + inv_S[0] * A_u[cell_id][0] * m_lambda_vector[i][0] +
+                          inv_S[1] * A_u[cell_id][1] * m_lambda_vector[i][1];
+        }
+      });
+    return M;
+  }
+
+  NodeArray<double>
+  compute_Flux1(const DiscreteFunctionP0Vector<const double>& F)
+  {
+    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const size_t nb_waves           = m_lambda_vector.size();
+
+    NodeArray<double> Flux_F{m_mesh.connectivity(), nb_waves};
+    Flux_F.fill(0);
+
+    parallel_for(
+      m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        for (size_t i = 0; i < nb_waves; ++i) {
+          const size_t N = m_mesh.numberOfCells();
+
+          const auto& cell_to_node   = cell_to_node_matrix[cell_id];
+          const NodeId left_node_id  = cell_to_node[0];
+          const NodeId right_node_id = cell_to_node[1];
+
+          const CellId prev_cell_id = (cell_id + N - 1) % N;
+          const CellId next_cell_id = (cell_id + 1) % N;
+
+          if (m_lambda_vector[i] < 0) {
+            Flux_F[left_node_id][i]  = F[cell_id][i];
+            Flux_F[right_node_id][i] = F[next_cell_id][i];
+
+          } else if (m_lambda_vector[i] > 0) {
+            Flux_F[left_node_id][i] = F[prev_cell_id][i];
+
+            Flux_F[right_node_id][i] = F[cell_id][i];
+
+          } else {
+            Flux_F[right_node_id][i] = 0;
+            Flux_F[left_node_id][i]  = 0;
+          }
+        }
+      });
+    return Flux_F;
+  }
+
+  template <size_t nb_subtimesteps>
+  std::vector<CellValue<double>>
+  compute_PFnp1(const DiscreteFunctionP0Vector<const double>& F0,
+                std::vector<DiscreteFunctionP0Vector<double>>& deltaF,
+                DiscreteFunctionP0Vector<double>& deltaF0,
+                const TinyMatrix<nb_subtimesteps>& B,
+                const TinyVector<nb_subtimesteps>& b0)
+
+  {
+    std::vector<CellValue<double>> PFnp1;
+    for (size_t i_subtimestep = 0; i_subtimestep < nb_subtimesteps; ++i_subtimestep) {
+      PFnp1.emplace_back(m_mesh.connectivity());
+    }
+
+    for (size_t i_subtimestep = 0; i_subtimestep < nb_subtimesteps; ++i_subtimestep) {
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+          const double inv_dx = m_inv_dx_table[cell_id];
+          const double c      = m_dt * inv_dx;
+
+          double test_rho_j = 0;
+
+          PFnp1[i_subtimestep][cell_id] = 0;
+          for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
+            const double lambda = m_lambda_vector[i_wave];
+            TinyVector<nb_subtimesteps> deltaF_vec;
+            TinyVector<nb_subtimesteps> b0_deltaF0;
+
+            for (size_t j_subtimestep = 0; j_subtimestep < nb_subtimesteps; ++j_subtimestep) {
+              deltaF_vec[j_subtimestep] = deltaF[j_subtimestep][cell_id][i_wave];
+              b0_deltaF0[j_subtimestep] = b0[j_subtimestep] * deltaF0[cell_id][i_wave];
+            }
+            const TinyVector<nb_subtimesteps> B_deltaF = B * deltaF_vec;
+
+            PFnp1[i_subtimestep][cell_id] +=
+              F0[cell_id][i_wave] - c * lambda * (B_deltaF[i_subtimestep] + b0_deltaF0[i_subtimestep]);
+
+            test_rho_j += -c * lambda * B_deltaF[i_subtimestep];
+          }
+        });
+    }
+    return PFnp1;
+  }
+
+  CellValue<double>
+  compute_PFn(const DiscreteFunctionP0Vector<const double> F)
+  {
+    CellValue<double> PFn{m_mesh.connectivity()};
+
+    // parallel_for(
+    //   m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+    for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+      double PF = 0;
+      for (size_t i = 0; i < m_lambda_vector.size(); ++i) {
+        PF += F[cell_id][i];
+      }
+      PFn[cell_id] = PF;
+    }   //    });
+    return PFn;
+  }
+
+  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply()
+  {
+    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    // const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+    const size_t nb_waves     = m_lambda_vector.size();
+    const size_t nb_equations = 4;
+
+    DiscreteFunctionP0Vector<double> F0_rho   = copy(m_Fn_rho);
+    DiscreteFunctionP0Vector<double> F0_rho_u = copy(m_Fn_rho_u1);
+    DiscreteFunctionP0Vector<double> F0_rho_E = copy(m_Fn_rho_E);
+
+    DiscreteFunctionP0Vector<double> F_rho(p_mesh, nb_waves);
+    DiscreteFunctionP0Vector<double> F_rho_u1(p_mesh, nb_waves);
+    DiscreteFunctionP0Vector<double> F_rho_u2(p_mesh, nb_waves);
+    DiscreteFunctionP0Vector<double> F_rho_E(p_mesh, nb_waves);
+
+    DiscreteFunctionP0Vector<double> Fnp1_rho(p_mesh, nb_waves);
+    DiscreteFunctionP0Vector<double> Fnp1_rho_u1(p_mesh, nb_waves);
+    DiscreteFunctionP0Vector<double> Fnp1_rho_u2(p_mesh, nb_waves);
+    DiscreteFunctionP0Vector<double> Fnp1_rho_E(p_mesh, nb_waves);
+
+    // Compute first order F
+
+    const CellValue<const double> rho   = compute_PFn(F0_rho);
+    const CellValue<const double> rho_u = compute_PFn(F0_rho_u);
+    const CellValue<const double> rho_E = compute_PFn(F0_rho_E);
+
+    return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(Fnp1_rho),
+                           std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u1),
+                           std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u2),
+                           std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_E));
+  }
+
+  EulerKineticSolver2D(std::shared_ptr<const MeshType> mesh,
+                       const double& dt,
+                       const double& eps,
+                       const double& gamma,
+                       const std::vector<TinyVector<MeshType::Dimention>>& lambda_vector,
+                       const size_t& SpaceOrder,
+                       const DiscreteFunctionP0Vector<const double>& Fn_rho,
+                       const DiscreteFunctionP0Vector<const double>& Fn_rho_u1,
+                       const DiscreteFunctionP0Vector<const double>& Fn_rho_u2,
+                       const DiscreteFunctionP0Vector<const double>& Fn_rho_E)
+    : m_dt{dt},
+      m_eps{eps},
+      m_lambda_vector{lambda_vector},
+      m_SpaceOrder{SpaceOrder},
+      m_Fn_rho{Fn_rho},
+      m_Fn_rho_u1{Fn_rho_u1},
+      m_Fn_rho_u2{Fn_rho_u2},
+      m_Fn_rho_E{Fn_rho_E},
+      m_gamma{gamma},
+      p_mesh{mesh},
+      m_mesh{*mesh}
+  {
+    if ((lambda_vector.size() != Fn_rho.size()) or (lambda_vector.size() != Fn_rho_u1.size()) or
+        (lambda_vector.size() != Fn_rho_u2.size()) or ((lambda_vector.size() != Fn_rho_E.size()))) {
+      throw NormalError("Incompatible dimensions of lambda vector and Fn");
+    }
+
+    m_inv_dx_table = [=] {
+      CellValue<double> inv_dx_table{m_mesh.connectivity()};
+      parallel_for(
+        m_mesh.numberOfCells(),
+        PUGS_LAMBDA(const CellId cell_id) { inv_dx_table[cell_id] = 1. / m_dx_table[cell_id]; });
+      return inv_dx_table;
+    }();
+  }
+
+  ~EulerKineticSolver2D() = default;
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolver2D(const double& dt,
+                     const double& gamma,
+                     const std::vector<TinyVector<MeshType::Dimention>>& lambda_vector,
+                     const double& eps,
+                     const size_t& SpaceOrder,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u1,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u2,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E)
+{
+  const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({F_rho, F_rho_u1, F_rho_u2, F_rho_E});
+  if (mesh_v.use_count() == 0) {
+    throw NormalError("F_rho, F_rho_u1, F_rho_u2 and F_rho_E have to be defined on the same mesh");
+  }
+  if (SpaceOrder > 1) {
+    throw NotImplementedError("Euler kinetic solver in 2D not implemented at order > 1");
+  }
+
+  return std::visit(
+    [&](auto&& p_mesh)
+      -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> {
+      using MeshType = mesh_type_t<decltype(p_mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 2)) {
+        auto Fn_rho    = F_rho->get<DiscreteFunctionP0Vector<const double>>();
+        auto Fn_rho_u1 = F_rho_u1->get<DiscreteFunctionP0Vector<const double>>();
+        auto Fn_rho_u2 = F_rho_u2->get<DiscreteFunctionP0Vector<const double>>();
+        auto Fn_rho_E  = F_rho_E->get<DiscreteFunctionP0Vector<const double>>();
+
+        EulerKineticSolver2D solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho, Fn_rho_u1, Fn_rho_u2,
+                                    Fn_rho_E);
+
+        return solver.apply();
+      } else {
+        throw NotImplementedError("Invalid mesh type for 2D Euler Kinetic solver");
+      }
+    },
+    mesh_v->variant());
+}

src/scheme/EulerKineticSolver2D.hpp

+#ifndef EULER_KINETIC_SOLVER_2D_HPP
+#define EULER_KINETIC_SOLVER_2D_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+getEulerKineticWaves2D(const std::vector<double>& lambda_vector,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho_u,
+                       const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
+                       const double& gamma);
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolver2D(const double& dt,
+                     const double& gamma,
+                     const std::vector<const TinyVector<MeshType::Dimention>>& lambda_vector,
+                     const double& eps,
+                     const size_t& SpaceOrder,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u1,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u2,
+                     const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E);
+
+#endif   // EULER_KINETIC_SOLVER_2D_HPP