From ef58044168d15a46ea694ea50708e32564533ccf Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Fri, 17 Jan 2025 11:09:37 +0100
Subject: [PATCH] Create kinetic scheme module

---
 src/language/modules/CMakeLists.txt          |   1 +
 src/language/modules/KineticSchemeModule.cpp | 929 +++++++++++++++++++
 src/language/modules/KineticSchemeModule.hpp |  23 +
 src/language/modules/SchemeModule.cpp        | 907 ------------------
 src/scheme/EulerKineticSolverMultiD.cpp      |  98 ++
 5 files changed, 1051 insertions(+), 907 deletions(-)
 create mode 100644 src/language/modules/KineticSchemeModule.cpp
 create mode 100644 src/language/modules/KineticSchemeModule.hpp

diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index b94d4249c..657693c1a 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -6,6 +6,7 @@ add_library(
   BuiltinModule.cpp
   CoreModule.cpp
   DevUtilsModule.cpp
+  KineticSchemeModule.cpp
   LinearSolverModule.cpp
   MathFunctionRegisterForVh.cpp
   MathModule.cpp
diff --git a/src/language/modules/KineticSchemeModule.cpp b/src/language/modules/KineticSchemeModule.cpp
new file mode 100644
index 000000000..a111e0522
--- /dev/null
+++ b/src/language/modules/KineticSchemeModule.cpp
@@ -0,0 +1,929 @@
+#include <language/modules/KineticSchemeModule.hpp>
+
+#include <language/modules/MeshModuleTypes.hpp>
+#include <language/modules/SchemeModuleTypes.hpp>
+
+#include <language/modules/MeshModuleTypes.hpp>
+#include <scheme/EulerKineticSolver.hpp>
+#include <scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.hpp>
+#include <scheme/EulerKineticSolverAcousticLagrangeFV.hpp>
+#include <scheme/EulerKineticSolverFirstOrderFV.hpp>
+#include <scheme/EulerKineticSolverLagrangeFV.hpp>
+#include <scheme/EulerKineticSolverMeanFluxMood.hpp>
+#include <scheme/EulerKineticSolverMoodFD.hpp>
+#include <scheme/EulerKineticSolverMoodFV.hpp>
+#include <scheme/EulerKineticSolverMultiD.hpp>
+#include <scheme/EulerKineticSolverMultiDOrder2.hpp>
+#include <scheme/EulerKineticSolverOneFluxMood.hpp>
+#include <scheme/EulerKineticSolverThirdOrderFD.hpp>
+#include <scheme/EulerKineticSolverThirdOrderFV.hpp>
+#include <scheme/EulerKineticSolverThirdOrderMoodFD.hpp>
+#include <scheme/EulerKineticSolverThirdOrderMoodFV.hpp>
+#include <scheme/Scalar2WavesKineticSolver.hpp>
+#include <scheme/Scalar3WavesKineticSolver.hpp>
+#include <scheme/ScalarKineticSolver.hpp>
+#include <scheme/WavesEquationKineticSolver.hpp>
+#include <scheme/WavesEquationKineticSolverThirdOrderTimeFD.hpp>
+#include <scheme/WavesEquationKineticSolverThirdOrderTimeFV.hpp>
+#include <scheme/WavesEquationNonUniformCellKineticSolver.hpp>
+#include <scheme/WavesEquationNonUniformKineticSolver.hpp>
+
+#include <language/modules/ModuleRepository.hpp>
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
+
+#include <memory>
+
+KineticSchemeModule::KineticSchemeModule()
+{
+  this->_addBuiltinFunction("scalar_2waves_kinetic",
+                            std::function(
+
+                              [](const double& dt, const FunctionSymbolId& flux_function_id, const double& a,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F)
+                                -> std::shared_ptr<const DiscreteFunctionVariant> {
+                                return scalar2WavesKineticSolver(dt, flux_function_id, a, eps, F);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("scalar_3waves_kinetic",
+                            std::function(
+
+                              [](const double& dt, const FunctionSymbolId& flux_function_id, const double& a,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F)
+                                -> std::shared_ptr<const DiscreteFunctionVariant> {
+                                return scalar3WavesKineticSolver(dt, flux_function_id, a, eps, F);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("scalar_kinetic", std::function(
+
+                                                [](const double& dt, const FunctionSymbolId& flux_function_id,
+                                                   const std::vector<double>& lambda_vector, const double& eps,
+                                                   const std::shared_ptr<const DiscreteFunctionVariant>& F)
+                                                  -> std::shared_ptr<const DiscreteFunctionVariant> {
+                                                  return scalarKineticSolver(dt, flux_function_id, lambda_vector, eps,
+                                                                             F);
+                                                }
+
+                                                ));
+
+  this->_addBuiltinFunction("get_scalar_kinetic_waves", std::function(
+
+                                                          [](const std::vector<double>& lambda_vector,
+                                                             const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                                             const std::shared_ptr<const DiscreteFunctionVariant>& Au)
+                                                            -> std::shared_ptr<const DiscreteFunctionVariant> {
+                                                            return getScalarKineticWaves(lambda_vector, u, Au);
+                                                          }));
+  this->_addBuiltinFunction("waves_equation_kinetic",
+                            std::function(
+
+                              [](const double& dt, const TinyMatrix<2>& A, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
+                                return wavesEquationKineticSolver(dt, A, lambda_vector, eps, F_u, F_p);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_waves_equation_kinetic_waves",
+                            std::function(
+
+                              [](const std::vector<double>& lambda_vector,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const TinyMatrix<2>& A) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                       std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getWavesEquationKineticWaves(lambda_vector, u, p, A);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("waves_equation_kinetic_third_order_time_FV",
+                            std::function(
+
+                              [](const double& dt, const TinyMatrix<2>& A, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
+                                return wavesEquationKineticSolverThirdOrderTimeFV(dt, A, lambda_vector, eps, F_u, F_p);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_waves_equation_kinetic_waves_third_order_time_FV",
+                            std::function(
+
+                              [](const std::vector<double>& lambda_vector,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const TinyMatrix<2>& A) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                       std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getWavesEquationKineticWavesThirdOrderTimeFV(lambda_vector, u, p, A);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("waves_equation_kinetic_third_order_time_FD",
+                            std::function(
+
+                              [](const double& dt, const TinyMatrix<2>& A, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
+                                return wavesEquationKineticSolverThirdOrderTimeFD(dt, A, lambda_vector, eps, F_u, F_p);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_waves_equation_kinetic_waves_third_order_time_FD",
+                            std::function(
+
+                              [](const std::vector<double>& lambda_vector,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
+                                 const TinyMatrix<2>& A) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                       std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getWavesEquationKineticWavesThirdOrderTimeFV(lambda_vector, u, p, A);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("waves_equation_non_uniform_kinetic",
+                            std::function(
+
+                              [](const double& dt, const TinyMatrix<2>& A,
+                                 const std::shared_ptr<const ItemArrayVariant>& lambda_vector, const double& eps,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
+                                return wavesEquationNonUniformKineticSolver(dt, A, lambda_vector, eps, F_u, F_p);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_waves_equation_non_uniform_kinetic_waves",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p, const TinyMatrix<2>& A,
+                                 const uint64_t& nb_waves,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const ItemArrayVariant>> {
+                                return getWavesEquationNonUniformKineticWaves(u, p, A, nb_waves, dt);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("waves_equation_non_uniform_cell_kinetic",
+                            std::function(
+
+                              [](const double& dt, const TinyMatrix<2>& A,
+                                 const std::shared_ptr<const ItemArrayVariant>& lambda_vector, const double& eps,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
+                                return wavesEquationNonUniformCellKineticSolver(dt, A, lambda_vector, eps, F_u, F_p);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_waves_equation_non_uniform_cell_kinetic_waves",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p, const TinyMatrix<2>& A,
+                                 const uint64_t& nb_waves,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const ItemArrayVariant>> {
+                                return getWavesEquationNonUniformCellKineticWaves(u, p, A, nb_waves, dt);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolver(dt, gamma, lambda_vector, eps, F_rho, F_rho_u, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWaves(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolution(t, rho_ex, u_ex);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_mean_flux_mood",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverMeanFluxMood(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
+                                                                      F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_mean_flux_mood",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesMeanFluxMood(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_mean_flux_mood",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolutionMeanFluxMood(t, rho_ex, u_ex);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_one_flux_mood",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverOneFluxMood(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
+                                                                     F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_one_flux_mood",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesOneFluxMood(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_one_flux_mood",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolutionOneFluxMood(t, rho_ex, u_ex);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_mood_FV",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverMoodFV(dt, gamma, lambda_vector, eps, F_rho, F_rho_u, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_mood_FV",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesMoodFV(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_mood_FV",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolutionMoodFV(t, rho_ex, u_ex);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_mood_FD",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverMoodFD(dt, gamma, lambda_vector, eps, F_rho, F_rho_u, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_mood_FD",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesMoodFD(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_mood_FD",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolutionMoodFD(t, rho_ex, u_ex);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_third_order_mood_FV",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverThirdOrderMoodFV(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
+                                                                          F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_mood_FV",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesThirdOrderMoodFV(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_mood_FV",
+                            std::function([](const double& t, const double& gamma,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return getEulerKineticExactSolutionThirdOrderMoodFV(t, gamma, rho_ex, u_ex, p_ex,
+                                                                                  rho_u_ex, rho_E_ex);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("euler_kinetic_third_order_mood_FD",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverThirdOrderMoodFD(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
+                                                                          F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_mood_FD",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesThirdOrderMoodFD(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_mood_FD",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolutionThirdOrderMoodFD(t, rho_ex, u_ex);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_third_order_FV",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const size_t& SpaceOrder, const size_t& Limitation,
+                                 const bool& MoodLimitation,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverThirdOrderFV(dt, gamma, lambda_vector, eps, SpaceOrder,
+                                                                      Limitation, MoodLimitation, F_rho, F_rho_u, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_FV",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesThirdOrderFV(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_FV",
+                            std::function([](const double& t, const double& gamma,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return getEulerKineticExactSolutionThirdOrderFV(t, gamma, rho_ex, u_ex, p_ex, rho_u_ex,
+                                                                              rho_E_ex);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("euler_kinetic_MultiD",
+                            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_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
+                                                                F_rho_u, F_E, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_MultiD",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& 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_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
+                                                                F_rho_u, F_E, bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_lambda_vector_2d",
+                            std::function(
+
+                              [](const std::shared_ptr<const MeshVariant>& mesh, const double& lambda, const size_t& L,
+                                 const size_t& M) -> std::vector<TinyVector<2>> {
+                                return getLambdaVector2D(mesh, lambda, L, M);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_lambda_vector_3d", std::function(
+
+                                                      [](const std::shared_ptr<const MeshVariant>& mesh,
+                                                         const double& lambda) -> std::vector<TinyVector<3>> {
+                                                        return getLambdaVector3D(mesh, lambda);
+                                                      }
+
+                                                      ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_MultiD",
+                            std::function(
+
+                              [](const std::vector<TinyVector<2>>& 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 getEulerKineticWavesMultiD(lambda_vector, rho, rho_u, rho_E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_MultiD",
+                            std::function(
+
+                              [](const std::vector<TinyVector<3>>& 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 getEulerKineticWavesMultiD(lambda_vector, rho, rho_u, rho_E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_MultiD_Order2",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector,
+                                 const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder,
+                                                                      TimeOrder, F_rho, F_rho_u, F_E,
+                                                                      bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_MultiD_Order2",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& lambda_vector,
+                                 const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder,
+                                                                      TimeOrder, F_rho, F_rho_u, F_E,
+                                                                      bc_descriptor_list);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("euler_kinetic_first_order_FV",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const size_t& SpaceOrder, const size_t& Limitation,
+                                 const bool& MoodLimitation,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverFirstOrderFV(dt, gamma, lambda_vector, eps, SpaceOrder,
+                                                                      Limitation, MoodLimitation, F_rho, F_rho_u, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_first_order_FV",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesFirstOrderFV(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_first_order_FV",
+                            std::function([](const double& t, const double& gamma,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return getEulerKineticExactSolutionFirstOrderFV(t, gamma, rho_ex, u_ex, p_ex, rho_u_ex,
+                                                                              rho_E_ex);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("refine",
+                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& old_u_old_mesh,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& new_u_new_mesh)
+                                            -> std::shared_ptr<const DiscreteFunctionVariant> {
+                              return refine(old_u_old_mesh, new_u_new_mesh);
+                            }));
+
+  this->_addBuiltinFunction("coarsen",
+                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& u_old_mesh,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& new_u_new_mesh)
+                                            -> std::shared_ptr<const DiscreteFunctionVariant> {
+                              return coarsen(u_old_mesh, new_u_new_mesh);
+                            }));
+
+  this->_addBuiltinFunction("euler_kinetic_lagrange_FV",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const size_t& space_order, const size_t& time_order, const double& eps,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E)
+                                -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverLagrangeFV(dt, gamma, lambda_vector, space_order, time_order,
+                                                                    eps, F_rho, F_rho_u, F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_lagrange_FV",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesLagrangeFV(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_lagrange_FV",
+                            std::function([](const double& t, const double& gamma,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return getEulerKineticExactSolutionLagrangeFV(t, gamma, rho_ex, u_ex, p_ex, rho_u_ex,
+                                                                            rho_E_ex);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV",
+                            std::function(
+
+                              [](const double& dt, const std::vector<double>& lambda_vector, const double& gamma,
+                                 const double& eps, const size_t& space_order, const size_t& time_order,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& p)
+                                -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return eulerKineticSolverAcousticLagrangeFV(dt, lambda_vector, gamma, eps, space_order,
+                                                                            time_order, rho, u, E, p);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_acoustic_lagrange_FV",
+                            std::function([](const double& t, const double& gamma,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return getEulerKineticExactSolutionAcousticLagrangeFV(t, gamma, rho_ex, u_ex, p_ex,
+                                                                                    rho_u_ex, rho_E_ex);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("euler_kinetic_acoustic2_lagrange_FV_order2",
+                            std::function([](const double& dt, const std::vector<double>& lambda_vector,
+                                             const double& gamma, const double& eps, const size_t& space_order,
+                                             const size_t& time_order, const bool& TVD_limitation,
+                                             const bool& MOOD_limitation,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p)
+                                            -> std::tuple<std::shared_ptr<const MeshVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return eulerKineticSolverAcoustic2LagrangeFVOrder2(dt, lambda_vector, gamma, eps,
+                                                                                 space_order, time_order,
+                                                                                 TVD_limitation, MOOD_limitation, rho,
+                                                                                 u, E, p);
+                            }));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_acoustic2_lagrange_FV_order2",
+                            std::function([](const double& t, const double& gamma,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
+                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>,
+                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
+                              return getEulerKineticExactSolutionAcoustic2LagrangeFVOrder2(t, gamma, rho_ex, u_ex, p_ex,
+                                                                                           rho_u_ex, rho_E_ex);
+                            }
+
+                                          ));
+
+  this->_addBuiltinFunction("euler_kinetic_third_order_FD",
+                            std::function(
+
+                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
+                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+                                 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 eulerKineticSolverThirdOrderFD(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
+                                                                      F_E);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_FD",
+                            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>& E,
+                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticWavesThirdOrderFD(lambda_vector, rho, rho_u, E, gamma);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_FD",
+                            std::function(
+                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
+                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                              std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return getEulerKineticExactSolutionThirdOrderFD(t, rho_ex, u_ex);
+                              }
+
+                              ));
+}
+
+void
+KineticSchemeModule::registerOperators() const
+{}
+
+void
+KineticSchemeModule::registerCheckpointResume() const
+{}
+
+ModuleRepository::Subscribe<KineticSchemeModule> kinetic_scheme_module;
diff --git a/src/language/modules/KineticSchemeModule.hpp b/src/language/modules/KineticSchemeModule.hpp
new file mode 100644
index 000000000..7881448e7
--- /dev/null
+++ b/src/language/modules/KineticSchemeModule.hpp
@@ -0,0 +1,23 @@
+#ifndef KINETIC_SCHEME_MODULE_HPP
+#define KINETIC_SCHEME_MODULE_HPP
+
+#include <language/modules/BuiltinModule.hpp>
+
+class KineticSchemeModule : public BuiltinModule
+{
+ public:
+  std::string_view
+  name() const final
+  {
+    return "kinetic_scheme";
+  }
+
+  void registerOperators() const final;
+  void registerCheckpointResume() const final;
+
+  KineticSchemeModule();
+
+  ~KineticSchemeModule() = default;
+};
+
+#endif   // KINETIC_SCHEME_MODULE_HPP
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 11cb53bc1..40cc703cc 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -18,8 +18,6 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/IBoundaryDescriptor.hpp>
 #include <mesh/IZoneDescriptor.hpp>
-#include <mesh/ItemArrayVariant.hpp>
-#include <mesh/ItemValueVariant.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -38,21 +36,6 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <scheme/DiscreteFunctionVectorIntegrator.hpp>
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
-#include <scheme/EulerKineticSolver.hpp>
-#include <scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.hpp>
-#include <scheme/EulerKineticSolverAcousticLagrangeFV.hpp>
-#include <scheme/EulerKineticSolverFirstOrderFV.hpp>
-#include <scheme/EulerKineticSolverLagrangeFV.hpp>
-#include <scheme/EulerKineticSolverMeanFluxMood.hpp>
-#include <scheme/EulerKineticSolverMoodFD.hpp>
-#include <scheme/EulerKineticSolverMoodFV.hpp>
-#include <scheme/EulerKineticSolverMultiD.hpp>
-#include <scheme/EulerKineticSolverMultiDOrder2.hpp>
-#include <scheme/EulerKineticSolverOneFluxMood.hpp>
-#include <scheme/EulerKineticSolverThirdOrderFD.hpp>
-#include <scheme/EulerKineticSolverThirdOrderFV.hpp>
-#include <scheme/EulerKineticSolverThirdOrderMoodFD.hpp>
-#include <scheme/EulerKineticSolverThirdOrderMoodFV.hpp>
 #include <scheme/ExternalBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/FluxingAdvectionSolver.hpp>
@@ -65,17 +48,9 @@
 #include <scheme/InflowListBoundaryConditionDescriptor.hpp>
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
-#include <scheme/Scalar2WavesKineticSolver.hpp>
-#include <scheme/Scalar3WavesKineticSolver.hpp>
-#include <scheme/ScalarKineticSolver.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
 #include <scheme/WallBoundaryConditionDescriptor.hpp>
-#include <scheme/WavesEquationKineticSolver.hpp>
-#include <scheme/WavesEquationKineticSolverThirdOrderTimeFD.hpp>
-#include <scheme/WavesEquationKineticSolverThirdOrderTimeFV.hpp>
-#include <scheme/WavesEquationNonUniformCellKineticSolver.hpp>
-#include <scheme/WavesEquationNonUniformKineticSolver.hpp>
 #include <utils/Socket.hpp>
 
 #include <utils/checkpointing/ReadDiscreteFunctionVariant.hpp>
@@ -435,888 +410,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("scalar_2waves_kinetic",
-                            std::function(
-
-                              [](const double& dt, const FunctionSymbolId& flux_function_id, const double& a,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F)
-                                -> std::shared_ptr<const DiscreteFunctionVariant> {
-                                return scalar2WavesKineticSolver(dt, flux_function_id, a, eps, F);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("scalar_3waves_kinetic",
-                            std::function(
-
-                              [](const double& dt, const FunctionSymbolId& flux_function_id, const double& a,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F)
-                                -> std::shared_ptr<const DiscreteFunctionVariant> {
-                                return scalar3WavesKineticSolver(dt, flux_function_id, a, eps, F);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("scalar_kinetic", std::function(
-
-                                                [](const double& dt, const FunctionSymbolId& flux_function_id,
-                                                   const std::vector<double>& lambda_vector, const double& eps,
-                                                   const std::shared_ptr<const DiscreteFunctionVariant>& F)
-                                                  -> std::shared_ptr<const DiscreteFunctionVariant> {
-                                                  return scalarKineticSolver(dt, flux_function_id, lambda_vector, eps,
-                                                                             F);
-                                                }
-
-                                                ));
-
-  this->_addBuiltinFunction("get_scalar_kinetic_waves", std::function(
-
-                                                          [](const std::vector<double>& lambda_vector,
-                                                             const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                                             const std::shared_ptr<const DiscreteFunctionVariant>& Au)
-                                                            -> std::shared_ptr<const DiscreteFunctionVariant> {
-                                                            return getScalarKineticWaves(lambda_vector, u, Au);
-                                                          }));
-  this->_addBuiltinFunction("waves_equation_kinetic",
-                            std::function(
-
-                              [](const double& dt, const TinyMatrix<2>& A, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
-                                return wavesEquationKineticSolver(dt, A, lambda_vector, eps, F_u, F_p);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_waves_equation_kinetic_waves",
-                            std::function(
-
-                              [](const std::vector<double>& lambda_vector,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
-                                 const TinyMatrix<2>& A) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                       std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getWavesEquationKineticWaves(lambda_vector, u, p, A);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("waves_equation_kinetic_third_order_time_FV",
-                            std::function(
-
-                              [](const double& dt, const TinyMatrix<2>& A, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
-                                return wavesEquationKineticSolverThirdOrderTimeFV(dt, A, lambda_vector, eps, F_u, F_p);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_waves_equation_kinetic_waves_third_order_time_FV",
-                            std::function(
-
-                              [](const std::vector<double>& lambda_vector,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
-                                 const TinyMatrix<2>& A) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                       std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getWavesEquationKineticWavesThirdOrderTimeFV(lambda_vector, u, p, A);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("waves_equation_kinetic_third_order_time_FD",
-                            std::function(
-
-                              [](const double& dt, const TinyMatrix<2>& A, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
-                                return wavesEquationKineticSolverThirdOrderTimeFD(dt, A, lambda_vector, eps, F_u, F_p);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_waves_equation_kinetic_waves_third_order_time_FD",
-                            std::function(
-
-                              [](const std::vector<double>& lambda_vector,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p,
-                                 const TinyMatrix<2>& A) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                       std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getWavesEquationKineticWavesThirdOrderTimeFV(lambda_vector, u, p, A);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("waves_equation_non_uniform_kinetic",
-                            std::function(
-
-                              [](const double& dt, const TinyMatrix<2>& A,
-                                 const std::shared_ptr<const ItemArrayVariant>& lambda_vector, const double& eps,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
-                                return wavesEquationNonUniformKineticSolver(dt, A, lambda_vector, eps, F_u, F_p);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_waves_equation_non_uniform_kinetic_waves",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p, const TinyMatrix<2>& A,
-                                 const uint64_t& nb_waves,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const ItemArrayVariant>> {
-                                return getWavesEquationNonUniformKineticWaves(u, p, A, nb_waves, dt);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("waves_equation_non_uniform_cell_kinetic",
-                            std::function(
-
-                              [](const double& dt, const TinyMatrix<2>& A,
-                                 const std::shared_ptr<const ItemArrayVariant>& lambda_vector, const double& eps,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_p)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>, double> {
-                                return wavesEquationNonUniformCellKineticSolver(dt, A, lambda_vector, eps, F_u, F_p);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_waves_equation_non_uniform_cell_kinetic_waves",
-                            std::function(
-
-                              [](const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p, const TinyMatrix<2>& A,
-                                 const uint64_t& nb_waves,
-                                 const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                 std::shared_ptr<const ItemArrayVariant>> {
-                                return getWavesEquationNonUniformCellKineticWaves(u, p, A, nb_waves, dt);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolver(dt, gamma, lambda_vector, eps, F_rho, F_rho_u, F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWaves(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolution(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_mean_flux_mood",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverMeanFluxMood(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
-                                                                      F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_mean_flux_mood",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesMeanFluxMood(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_mean_flux_mood",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolutionMeanFluxMood(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_one_flux_mood",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverOneFluxMood(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
-                                                                     F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_one_flux_mood",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesOneFluxMood(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_one_flux_mood",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolutionOneFluxMood(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_mood_FV",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverMoodFV(dt, gamma, lambda_vector, eps, F_rho, F_rho_u, F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_mood_FV",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesMoodFV(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_mood_FV",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolutionMoodFV(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_mood_FD",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverMoodFD(dt, gamma, lambda_vector, eps, F_rho, F_rho_u, F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_mood_FD",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesMoodFD(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_mood_FD",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolutionMoodFD(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_third_order_mood_FV",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverThirdOrderMoodFV(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
-                                                                          F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_mood_FV",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesThirdOrderMoodFV(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_mood_FV",
-                            std::function([](const double& t, const double& gamma,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
-                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return getEulerKineticExactSolutionThirdOrderMoodFV(t, gamma, rho_ex, u_ex, p_ex,
-                                                                                  rho_u_ex, rho_E_ex);
-                            }
-
-                                          ));
-
-  this->_addBuiltinFunction("euler_kinetic_third_order_mood_FD",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverThirdOrderMoodFD(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
-                                                                          F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_mood_FD",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesThirdOrderMoodFD(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_mood_FD",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolutionThirdOrderMoodFD(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_third_order_FV",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const size_t& SpaceOrder, const size_t& Limitation,
-                                 const bool& MoodLimitation,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverThirdOrderFV(dt, gamma, lambda_vector, eps, SpaceOrder,
-                                                                      Limitation, MoodLimitation, F_rho, F_rho_u, F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_FV",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesThirdOrderFV(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_FV",
-                            std::function([](const double& t, const double& gamma,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
-                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return getEulerKineticExactSolutionThirdOrderFV(t, gamma, rho_ex, u_ex, p_ex, rho_u_ex,
-                                                                              rho_E_ex);
-                            }
-
-                                          ));
-
-  this->_addBuiltinFunction("euler_kinetic_MultiD",
-                            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_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
-                                                                F_rho_u, F_E, bc_descriptor_list);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_MultiD",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& 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_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
-                                                                F_rho_u, F_E, bc_descriptor_list);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_lambda_vector_2d",
-                            std::function(
-
-                              [](const std::shared_ptr<const MeshVariant>& mesh, const double& lambda, const size_t& L,
-                                 const size_t& M) -> std::vector<TinyVector<2>> {
-                                return getLambdaVector2D(mesh, lambda, L, M);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_lambda_vector_3d", std::function(
-
-                                                      [](const std::shared_ptr<const MeshVariant>& mesh,
-                                                         const double& lambda) -> std::vector<TinyVector<3>> {
-                                                        return getLambdaVector3D(mesh, lambda);
-                                                      }
-
-                                                      ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_MultiD",
-                            std::function(
-
-                              [](const std::vector<TinyVector<2>>& 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 getEulerKineticWavesMultiD(lambda_vector, rho, rho_u, rho_E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_MultiD",
-                            std::function(
-
-                              [](const std::vector<TinyVector<3>>& 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 getEulerKineticWavesMultiD(lambda_vector, rho, rho_u, rho_E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_MultiD_Order2",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector,
-                                 const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder,
-                                                                      TimeOrder, F_rho, F_rho_u, F_E,
-                                                                      bc_descriptor_list);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_MultiD_Order2",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& lambda_vector,
-                                 const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                     std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder,
-                                                                      TimeOrder, F_rho, F_rho_u, F_E,
-                                                                      bc_descriptor_list);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("euler_kinetic_first_order_FV",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const size_t& SpaceOrder, const size_t& Limitation,
-                                 const bool& MoodLimitation,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverFirstOrderFV(dt, gamma, lambda_vector, eps, SpaceOrder,
-                                                                      Limitation, MoodLimitation, F_rho, F_rho_u, F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_first_order_FV",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesFirstOrderFV(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_first_order_FV",
-                            std::function([](const double& t, const double& gamma,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
-                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return getEulerKineticExactSolutionFirstOrderFV(t, gamma, rho_ex, u_ex, p_ex, rho_u_ex,
-                                                                              rho_E_ex);
-                            }
-
-                                          ));
-
-  this->_addBuiltinFunction("refine",
-                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& old_u_old_mesh,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& new_u_new_mesh)
-                                            -> std::shared_ptr<const DiscreteFunctionVariant> {
-                              return refine(old_u_old_mesh, new_u_new_mesh);
-                            }));
-
-  this->_addBuiltinFunction("coarsen",
-                            std::function([](const std::shared_ptr<const DiscreteFunctionVariant>& u_old_mesh,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& new_u_new_mesh)
-                                            -> std::shared_ptr<const DiscreteFunctionVariant> {
-                              return coarsen(u_old_mesh, new_u_new_mesh);
-                            }));
-
-  this->_addBuiltinFunction("euler_kinetic_lagrange_FV",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const size_t& space_order, const size_t& time_order, const double& eps,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_E)
-                                -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverLagrangeFV(dt, gamma, lambda_vector, space_order, time_order,
-                                                                    eps, F_rho, F_rho_u, F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_lagrange_FV",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesLagrangeFV(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_lagrange_FV",
-                            std::function([](const double& t, const double& gamma,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
-                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return getEulerKineticExactSolutionLagrangeFV(t, gamma, rho_ex, u_ex, p_ex, rho_u_ex,
-                                                                            rho_E_ex);
-                            }
-
-                                          ));
-
-  this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV",
-                            std::function(
-
-                              [](const double& dt, const std::vector<double>& lambda_vector, const double& gamma,
-                                 const double& eps, const size_t& space_order, const size_t& time_order,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& p)
-                                -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return eulerKineticSolverAcousticLagrangeFV(dt, lambda_vector, gamma, eps, space_order,
-                                                                            time_order, rho, u, E, p);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_acoustic_lagrange_FV",
-                            std::function([](const double& t, const double& gamma,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
-                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return getEulerKineticExactSolutionAcousticLagrangeFV(t, gamma, rho_ex, u_ex, p_ex,
-                                                                                    rho_u_ex, rho_E_ex);
-                            }
-
-                                          ));
-
-  this->_addBuiltinFunction("euler_kinetic_acoustic2_lagrange_FV_order2",
-                            std::function([](const double& dt, const std::vector<double>& lambda_vector,
-                                             const double& gamma, const double& eps, const size_t& space_order,
-                                             const size_t& time_order, const bool& TVD_limitation,
-                                             const bool& MOOD_limitation,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& E,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p)
-                                            -> std::tuple<std::shared_ptr<const MeshVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return eulerKineticSolverAcoustic2LagrangeFVOrder2(dt, lambda_vector, gamma, eps,
-                                                                                 space_order, time_order,
-                                                                                 TVD_limitation, MOOD_limitation, rho,
-                                                                                 u, E, p);
-                            }));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_acoustic2_lagrange_FV_order2",
-                            std::function([](const double& t, const double& gamma,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& p_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
-                                             const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
-                                            -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>,
-                                                          std::shared_ptr<const DiscreteFunctionVariant>> {
-                              return getEulerKineticExactSolutionAcoustic2LagrangeFVOrder2(t, gamma, rho_ex, u_ex, p_ex,
-                                                                                           rho_u_ex, rho_E_ex);
-                            }
-
-                                          ));
-
-  this->_addBuiltinFunction("euler_kinetic_third_order_FD",
-                            std::function(
-
-                              [](const double& dt, const double& gamma, const std::vector<double>& lambda_vector,
-                                 const double& eps, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
-                                 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 eulerKineticSolverThirdOrderFD(dt, gamma, lambda_vector, eps, F_rho, F_rho_u,
-                                                                      F_E);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_waves_third_order_FD",
-                            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>& E,
-                                 const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>,
-                                                                    std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticWavesThirdOrderFD(lambda_vector, rho, rho_u, E, gamma);
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("get_euler_kinetic_exact_solution_third_order_FD",
-                            std::function(
-                              [](const double& t, const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& u_ex)
-                                -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
-                                              std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return getEulerKineticExactSolutionThirdOrderFD(t, rho_ex, u_ex);
-                              }
-
-                              ));
-
   this->_addBuiltinFunction("inflow", std::function(
 
                                         [](std::shared_ptr<const IBoundaryDescriptor> boundary,
diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index bc2fe33ba..791cf9a06 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -1077,6 +1077,104 @@ class EulerKineticSolverMultiD
             NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
             nr.fill(zero);
 
+            for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+              const NodeId node_id                      = node_list[i_node];
+              const auto& node_to_cell                  = node_to_cell_matrix[node_id];
+              const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+
+              for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+                const CellId cell_id     = node_to_cell[i_cell];
+                const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
+
+                nr[node_id] += njr(cell_id, i_node_cell);
+              }
+              nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
+              auto nxn    = tensorProduct(nr[node_id], nr[node_id]);
+              auto txt    = I - nxn;
+
+              for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+                const auto& node_to_face                  = node_to_face_matrix[node_id];
+                const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
+
+                double sum                = 0;
+                Fr_rho[node_id][i_wave]   = 0;
+                Fr_rho_u[node_id][i_wave] = zero;
+                Fr_rho_E[node_id][i_wave] = 0;
+
+                for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
+                  const FaceId face_id                      = node_to_face[i_face];
+                  const auto& face_to_cell                  = face_to_cell_matrix[face_id];
+                  const size_t i_node_face                  = node_local_number_in_its_face[i_face];
+                  const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+
+                  for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+                    const CellId cell_id     = face_to_cell[i_cell];
+                    const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+                    const double sign            = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+                    TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face);
+
+                    double li_nlr = dot(m_lambda_vector[i_wave], n_face);
+
+                    if (li_nlr > 0) {
+                      Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
+                      Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
+
+                      sum += li_nlr;
+                    }
+
+                    TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+                    double li_nlr_prime             = dot(m_lambda_vector[i_wave], nlr_prime);
+
+                    if (li_nlr_prime > 0) {
+                      double rhoj_prime                  = rho[cell_id];
+                      TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
+                      double rho_Ej_prime                = rho_E[cell_id];
+
+                      double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
+                      double p     = (m_gamma - 1) * rho_e;
+
+                      TinyVector<Dimension> A_rho_prime = rho_uj_prime;
+                      TinyMatrix<Dimension> A_rho_u_prime =
+                        1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
+                      TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
+
+                      double Fn_rho_prime                  = rhoj_prime / nb_waves;
+                      TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
+                      double Fn_rho_E_prime                = rho_Ej_prime / nb_waves;
+
+                      for (size_t d1 = 0; d1 < Dimension; ++d1) {
+                        Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
+                        for (size_t d2 = 0; d2 < Dimension; ++d2) {
+                          Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
+                        }
+                        Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
+                      }
+
+                      Fr_rho[node_id][i_wave] += Fn_rho_prime * li_nlr_prime;
+                      Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u_prime;
+                      Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_nlr_prime;
+
+                      sum += li_nlr_prime;
+                    }
+                  }
+                }
+                if (sum != 0) {
+                  Fr_rho[node_id][i_wave] /= sum;
+                  Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+                  Fr_rho_E[node_id][i_wave] /= sum;
+                }
+              }
+            }
+          } else if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+            auto node_list = bc.nodeList();
+
+            TinyMatrix<Dimension> I = identity;
+
+            NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+            nr.fill(zero);
+
             for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
               const NodeId node_id                      = node_list[i_node];
               const auto& node_to_cell                  = node_to_cell_matrix[node_id];
-- 
GitLab