diff --git a/src/language/modules/KineticSchemeModule.cpp b/src/language/modules/KineticSchemeModule.cpp
index de1b5fd06d92c94ffa1b892c00b6a593bd45d36c..8a1f13b1e2c4d1b0afd2392a62b1687fc1ba4053 100644
--- a/src/language/modules/KineticSchemeModule.cpp
+++ b/src/language/modules/KineticSchemeModule.cpp
@@ -12,6 +12,7 @@
#include <scheme/EulerKineticSolverLagrangeMultiD.hpp>
#include <scheme/EulerKineticSolverLagrangeMultiDLocal.hpp>
#include <scheme/EulerKineticSolverLagrangeMultiDLocalOrder2.hpp>
+#include <scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.hpp>
#include <scheme/EulerKineticSolverMeanFluxMood.hpp>
#include <scheme/EulerKineticSolverMoodFD.hpp>
#include <scheme/EulerKineticSolverMoodFV.hpp>
@@ -904,6 +905,27 @@ KineticSchemeModule::KineticSchemeModule()
));
+ this->_addBuiltinFunction("euler_kinetic_lagrange_multiD_local_order2_periodic",
+ std::function(
+
+ [](const double& dt, const size_t& L, const size_t& M, const size_t& k,
+ 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,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return eulerKineticSolverLagrangeMultiDLocalOrder2(dt, L, M, k, gamma, eps, space_order,
+ time_order, rho, u, E, p,
+ bc_descriptor_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV",
std::function(
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index cea3033cd68b47769e56f7c97529aabe3b1fa471..c43dcecd3845c31d8d6ae757f22d819acdda3442 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -20,6 +20,7 @@ add_library(
EulerKineticSolverLagrangeMultiD.cpp
EulerKineticSolverLagrangeMultiDLocal.cpp
EulerKineticSolverLagrangeMultiDLocalOrder2.cpp
+ EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.cpp
EulerKineticSolverMeanFluxMood.cpp
EulerKineticSolverOneFluxMood.cpp
EulerKineticSolverMoodFD.cpp
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.cpp b/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..46682295729de63aa20b00505e66ad9737e28a25
--- /dev/null
+++ b/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.cpp
@@ -0,0 +1,666 @@
+#include <scheme/EulerKineticSolverLagrangeMultiDLocalOrder2.hpp>
+
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/StencilManager.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <tuple>
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
+{
+ private:
+ constexpr static size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+
+ const double m_dt;
+ const size_t m_L;
+ const size_t m_M;
+ const size_t m_k;
+ const double m_gamma;
+ const size_t m_spaceOrder;
+ const size_t m_timeOrder;
+ const CellValue<const double> m_rho;
+ const CellValue<const TinyVector<Dimension>> m_u;
+ const CellValue<const double> m_E;
+ const CellValue<const double> m_p;
+ const std::shared_ptr<const MeshType> p_mesh;
+ const MeshType& m_mesh;
+ const CellValue<const double> m_Vj = MeshDataManager::instance().getMeshData(m_mesh).Vj();
+ CellValue<const double> m_inv_Mj;
+ CellValue<const double> m_inv_Vj;
+ CellValue<const double> m_Mj;
+ const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+ const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+
+ public:
+ NodeArray<double>
+ compute_Flux_Fp(const DiscreteFunctionDPk<Dimension, const double> DPk_p,
+ const DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>> DPk_u,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const size_t nb_waves = m_k;
+ NodeArray<double> Fr(m_mesh.connectivity(), nb_waves);
+ const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+ const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
+
+ Fr.fill(0.);
+
+ parallel_for(
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ TinyVector<Dimension> inv_S = zero;
+ for (size_t d = 0; d < Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += lambda_vector[node_id][i][d] * lambda_vector[node_id][i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+ double lambda_r = std::sqrt(dot(lambda_vector[node_id][m_L], lambda_vector[node_id][m_L]));
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ double sum_li_njr = 0;
+
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+
+ double li_njr = dot(lambda_vector[node_id][i_wave], njr(cell_id, i_node));
+ if (li_njr > 1e-14) {
+ double Fj = (1. / nb_waves) * DPk_p[cell_id](m_xr[node_id]);
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fj += inv_S[d] * lambda_vector[node_id][i_wave][d] *
+ (lambda_r * lambda_r * DPk_u[cell_id](m_xr[node_id])[d]);
+ }
+
+ Fr[node_id][i_wave] += li_njr * Fj;
+ sum_li_njr += li_njr;
+ }
+ }
+ if (sum_li_njr != 0) {
+ Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
+ }
+ }
+ });
+
+ return Fr;
+ }
+
+ NodeArray<TinyVector<Dimension>>
+ compute_Flux_Fu(const DiscreteFunctionDPk<Dimension, const double> DPk_p,
+ const DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>> DPk_u,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const size_t nb_waves = m_k;
+ NodeArray<TinyVector<Dimension>> Fr(m_mesh.connectivity(), nb_waves);
+ const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+ const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
+
+ Fr.fill(zero);
+
+ parallel_for(
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ TinyVector<Dimension> inv_S = zero;
+ for (size_t d = 0; d < Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += lambda_vector[node_id][i][d] * lambda_vector[node_id][i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+
+ double sum_li_njr = 0;
+
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+ const unsigned int i_node = node_local_number_in_its_cell[i_cell];
+
+ double li_njr = dot(lambda_vector[node_id][i_wave], njr(cell_id, i_node));
+ if (li_njr > 1e-14) {
+ TinyVector<Dimension> Fj = (1. / nb_waves) * DPk_u[cell_id](m_xr[node_id]);
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fj[d] += 1. / Dimension * inv_S[d] * lambda_vector[node_id][i_wave][d] * DPk_p[cell_id](m_xr[node_id]);
+ }
+
+ Fr[node_id][i_wave] += li_njr * Fj;
+ sum_li_njr += li_njr;
+ }
+ }
+ if (sum_li_njr != 0) {
+ Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
+ }
+ }
+ });
+
+ return Fr;
+ }
+
+ const NodeValue<TinyVector<Dimension>>
+ compute_velocity_flux(NodeArray<double>& F, const NodeArray<const TinyVector<Dimension>>& lambda_vector)
+ {
+ NodeValue<TinyVector<Dimension>> u{p_mesh->connectivity()};
+ const size_t nb_waves = m_k;
+ u.fill(zero);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const double inv_lambda_squared = 1. / dot(lambda_vector[node_id][m_L], lambda_vector[node_id][m_L]);
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ u[node_id] += F[node_id][i_wave] * lambda_vector[node_id][i_wave];
+ }
+ u[node_id] = inv_lambda_squared * u[node_id];
+ });
+ return u;
+ }
+
+ const NodeValue<TinyVector<Dimension>>
+ compute_pressure_flux(NodeArray<TinyVector<Dimension>>& F,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector)
+ {
+ NodeValue<TinyVector<Dimension>> p{p_mesh->connectivity()};
+ const size_t nb_waves = m_k;
+ p.fill(zero);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ p[node_id][dim] += dot(F[node_id][i_wave], lambda_vector[node_id][i_wave]);
+ }
+ }
+ });
+ return p;
+ }
+
+ const CellValue<const double>
+ compute_delta_velocity(const NodeValue<const TinyVector<Dimension>>& ur)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+ CellValue<double> delta{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ delta.fill(0.);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+ for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id = cell_to_node[i_node];
+ delta[cell_id] += dot(Cjr(cell_id, i_node), ur[node_id]);
+ }
+ });
+ return delta;
+ }
+
+ const CellValue<const TinyVector<Dimension>>
+ compute_delta_pressure(const NodeValue<const TinyVector<Dimension>>& pr)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+ CellValue<TinyVector<Dimension>> delta{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ delta.fill(zero);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+ for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ const NodeId node_id = cell_to_node[i_node];
+ delta[cell_id][dim] += Cjr(cell_id, i_node)[dim] * pr[node_id][dim];
+ }
+ }
+ });
+ return delta;
+ }
+
+ const CellValue<const double>
+ compute_delta_pu(const NodeValue<const TinyVector<Dimension>>& ur, const NodeValue<const TinyVector<Dimension>>& pr)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+ CellValue<double> delta{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ delta.fill(0.);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+
+ for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id = cell_to_node[i_node];
+ TinyVector<Dimension> Fjr = zero;
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Fjr[dim] = Cjr(cell_id, i_node)[dim] * pr[node_id][dim];
+ }
+ delta[cell_id] += dot(Fjr, ur[node_id]);
+ }
+ });
+ return delta;
+ }
+
+ const CellValue<const double>
+ build_lambda_j(const CellValue<double>& rho, const CellValue<double>& p)
+ {
+ CellValue<double> lambda_j{p_mesh->connectivity()};
+ CellValue<double> cj{m_mesh.connectivity()};
+
+ if constexpr (Dimension == 1) {
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]);
+ lambda_j[cell_id] = std::sqrt((3. * (m_k - 1.)) / (m_k + 1.)) * std::abs(rho[cell_id] * cj[cell_id]);
+ });
+ } else if constexpr (Dimension == 2) {
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]);
+ lambda_j[cell_id] =
+ (12. * m_L * m_L) / ((m_L + 1.) * (2. * m_L + 1.)) * 2. * std::abs(rho[cell_id] * cj[cell_id]);
+ });
+ } else if constexpr (Dimension == 3) {
+ if (m_k != 6) {
+ throw NotImplementedError("3D Lagrangian scheme not implemented for other than 6 waves model");
+ }
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]);
+ lambda_j[cell_id] = 3. * std::abs(rho[cell_id] * cj[cell_id]);
+ });
+ }
+ return lambda_j;
+ }
+
+ const NodeArray<const TinyVector<Dimension>>
+ build_nodal_wavespeeds(const CellValue<const double>& lambda_j)
+ {
+ const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ NodeArray<TinyVector<Dimension>> lambda_vector(p_mesh->connectivity(), m_k);
+ const double pi = std::acos(-1);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ double lambda_r = 0.;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+
+ lambda_r = std::max(lambda_r, lambda_j[cell_id]);
+
+ // std::cout << "node_id = " << node_id << ", lambda_r = " << lambda_r << "\n";
+ }
+
+ if constexpr (Dimension == 1) {
+ for (size_t i_wave = 0; i_wave < m_k; ++i_wave) {
+ lambda_vector[node_id][i_wave][0] = lambda_r * (1 - 2. * i_wave / (m_k - 1));
+ }
+ } else if constexpr (Dimension == 2) {
+ for (size_t l = 1; l < m_L + 1; ++l) {
+ for (size_t m = 1; m < 4 * m_M + 1; ++m) {
+ size_t i_wave = (l - 1) * 4 * m_M + m - 1;
+ double ll = l;
+ lambda_vector[node_id][i_wave] =
+ TinyVector<Dimension>{(ll / m_L) * lambda_r * std::cos((m * pi) / (2 * m_M)), // + 0.2 * pi),
+ (ll / m_L) * lambda_r * std::sin((m * pi) / (2 * m_M))}; // + 0.2 * pi)};
+ }
+ }
+ } else if constexpr (Dimension == 3) {
+ if (m_k != 6) {
+ throw NotImplementedError("3D Lagrangian scheme not implemented for other than 6 waves model");
+ }
+ lambda_vector[node_id][0] = TinyVector<Dimension>{lambda_r, 0., 0.};
+ lambda_vector[node_id][1] = TinyVector<Dimension>{-lambda_r, 0., 0.};
+ lambda_vector[node_id][2] = TinyVector<Dimension>{0., lambda_r, 0.};
+ lambda_vector[node_id][3] = TinyVector<Dimension>{0., -lambda_r, 0.};
+ lambda_vector[node_id][4] = TinyVector<Dimension>{0., 0., lambda_r};
+ lambda_vector[node_id][5] = TinyVector<Dimension>{0., 0., -lambda_r};
+ }
+ });
+ return lambda_vector;
+ }
+
+ void
+ scalar_limiter(const CellValue<const double>& u, DiscreteFunctionDPk<Dimension, double>& DPk_uh) const
+ {
+ StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+ auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor);
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double uj = u[cell_id];
+
+ double u_min = uj;
+ double u_max = uj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ u_min = std::min(u_min, u[cell_stencil[i_cell]]);
+ u_max = std::max(u_max, u[cell_stencil[i_cell]]);
+ }
+
+ const double eps = 1E-14;
+ double coef = 1;
+ double lambda = 1.;
+
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ const CellId cell_i_id = cell_stencil[i_cell];
+ const double u_bar_i = DPk_uh[cell_id](m_xj[cell_i_id]);
+ const double Delta = (u_bar_i - uj);
+
+ if (Delta > eps) {
+ coef = std::min(1., (u_max - uj) / Delta);
+ } else if (Delta < -eps) {
+ coef = std::min(1., (u_min - uj) / Delta);
+ }
+ lambda = std::min(lambda, coef);
+ }
+
+ auto coefficients = DPk_uh.coefficients(cell_id);
+ coefficients[0] = (1 - lambda) * u[cell_id] + lambda * coefficients[0];
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda;
+ }
+ });
+ }
+
+ void
+ vector_limiter(const CellValue<const TinyVector<Dimension>>& u,
+ DiscreteFunctionDPk<Dimension, TinyVector<Dimension>>& DPk_uh) const
+ {
+ StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+ auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor);
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const TinyVector<Dimension> uj = u[cell_id];
+
+ TinyVector<Dimension> u_min = uj;
+ TinyVector<Dimension> u_max = uj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ u_min[dim] = std::min(u_min[dim], u[cell_stencil[i_cell]][dim]);
+ u_max[dim] = std::max(u_max[dim], u[cell_stencil[i_cell]][dim]);
+ }
+ }
+
+ const double eps = 1E-14;
+ TinyVector<Dimension> coef;
+ double lambda = 1.;
+
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ const CellId cell_i_id = cell_stencil[i_cell];
+ const TinyVector<Dimension> u_bar_i = DPk_uh[cell_id](m_xj[cell_i_id]);
+ TinyVector<Dimension> Delta;
+
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Delta[dim] = (u_bar_i[dim] - uj[dim]);
+ coef[dim] = 1.;
+ if (Delta[dim] > eps) {
+ coef[dim] = std::min(1., (u_max[dim] - uj[dim]) / Delta[dim]);
+ } else if (Delta[dim] < -eps) {
+ coef[dim] = std::min(1., (u_min[dim] - uj[dim]) / Delta[dim]);
+ }
+ lambda = std::min(lambda, coef[dim]);
+ }
+ }
+ auto coefficients = DPk_uh.coefficients(cell_id);
+
+ coefficients[0] = (1 - lambda) * u[cell_id] + lambda * coefficients[0];
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda;
+ }
+ });
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply()
+ {
+ const bool limitation = true;
+
+ const CellValue<double> rho = copy(m_rho);
+ const CellValue<TinyVector<Dimension>> u = copy(m_u);
+ const CellValue<double> p = copy(m_p);
+ const CellValue<double> E = copy(m_E);
+
+ CellValue<TinyVector<Dimension>> u_np1 = copy(u);
+ CellValue<double> p_np1 = copy(p);
+ CellValue<double> E_np1 = copy(E);
+
+ const CellValue<const double> lambda_j = build_lambda_j(rho, p);
+ const NodeArray<const TinyVector<Dimension>> lambda_vector = build_nodal_wavespeeds(lambda_j);
+
+ PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1);
+
+ auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, p),
+ DiscreteFunctionP0(p_mesh, u));
+
+ auto DPk_p = reconstruction[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto DPk_u = reconstruction[1]->template get<DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>>>();
+
+ DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_p);
+ DiscreteFunctionDPk<Dimension, TinyVector<Dimension>> u_lim = copy(DPk_u);
+
+ if (limitation) {
+ scalar_limiter(p, p_lim);
+ vector_limiter(u, u_lim);
+
+ synchronize(p_lim.cellArrays());
+ synchronize(u_lim.cellArrays());
+ }
+
+ NodeArray<double> Fr_p = compute_Flux_Fp(p_lim, u_lim, lambda_vector);
+ NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux_Fu(p_lim, u_lim, lambda_vector);
+
+ const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda_vector);
+ const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u, lambda_vector);
+
+ const CellValue<const double> delta_u = compute_delta_velocity(ur);
+ const CellValue<const TinyVector<Dimension>> delta_p = compute_delta_pressure(pr);
+ const CellValue<const double> delta_pu = compute_delta_pu(ur, pr);
+
+ CellValue<double> pj_star{m_mesh.connectivity()};
+ CellValue<TinyVector<Dimension>> uj_star = copy(m_u);
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double dt_over_Mj = m_dt * m_inv_Mj[cell_id];
+ double tauj_star = 1. / m_rho[cell_id] + dt_over_Mj * delta_u[cell_id];
+ pj_star[cell_id] = m_p[cell_id] / std::pow(m_rho[cell_id] * tauj_star, m_gamma);
+ uj_star[cell_id] += dt_over_Mj * delta_p[cell_id];
+ });
+
+ auto reconstruction_star =
+ PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, pj_star),
+ DiscreteFunctionP0(p_mesh, uj_star));
+
+ auto DPk_p_star = reconstruction_star[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto DPk_u_star =
+ reconstruction_star[1]->template get<DiscreteFunctionDPk<Dimension, const TinyVector<Dimension>>>();
+
+ DiscreteFunctionDPk<Dimension, double> p_star_lim = copy(DPk_p_star);
+ DiscreteFunctionDPk<Dimension, TinyVector<Dimension>> u_star_lim = copy(DPk_u_star);
+
+ if (limitation) {
+ scalar_limiter(pj_star, p_star_lim);
+ vector_limiter(uj_star, u_star_lim);
+
+ synchronize(p_star_lim.cellArrays());
+ synchronize(u_star_lim.cellArrays());
+ }
+
+ NodeArray<double> Fr_p_star = compute_Flux_Fp(p_star_lim, u_star_lim, lambda_vector);
+ NodeArray<TinyVector<Dimension>> Fr_u_star = compute_Flux_Fu(p_star_lim, u_star_lim, lambda_vector);
+
+ const NodeValue<TinyVector<Dimension>> ur_star = compute_velocity_flux(Fr_p_star, lambda_vector);
+ const NodeValue<TinyVector<Dimension>> pr_star = compute_pressure_flux(Fr_u_star, lambda_vector);
+
+ const CellValue<const double> delta_u_star = compute_delta_velocity(ur_star);
+ const CellValue<const TinyVector<Dimension>> delta_p_star = compute_delta_pressure(pr_star);
+ const CellValue<const double> delta_pu_star = compute_delta_pu(ur_star, pr_star);
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double dt_over_Mj = m_dt * m_inv_Mj[cell_id];
+ u_np1[cell_id] -= 0.5 * dt_over_Mj * (delta_p[cell_id] + delta_p_star[cell_id]);
+ E_np1[cell_id] -= 0.5 * dt_over_Mj * (delta_pu[cell_id] + delta_pu_star[cell_id]);
+ });
+
+ NodeValue<TinyVector<Dimension>> new_xr = copy(m_mesh.xr());
+ parallel_for(
+ m_mesh.numberOfNodes(),
+ PUGS_LAMBDA(NodeId node_id) { new_xr[node_id] += 0.5 * m_dt * (ur[node_id] + ur_star[node_id]); });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), new_xr);
+
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+ CellValue<double> rho_np1{p_mesh->connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rho_np1[cell_id] = m_Mj[cell_id] / new_Vj[cell_id]; });
+
+ return std::make_tuple(std::make_shared<MeshVariant>(new_mesh),
+ std::make_shared<DiscreteFunctionVariant>(
+ DiscreteFunctionP0<const double>(new_mesh, rho_np1)),
+ std::make_shared<DiscreteFunctionVariant>(
+ DiscreteFunctionP0<const TinyVector<Dimension>>(new_mesh, u_np1)),
+ std::make_shared<DiscreteFunctionVariant>(
+ DiscreteFunctionP0<const double>(new_mesh, E_np1)));
+ }
+
+ EulerKineticSolverLagrangeMultiDLocalOrder2Periodic(std::shared_ptr<const MeshType> mesh,
+ const double& dt,
+ const size_t& L,
+ const size_t& M,
+ const size_t& k,
+ const double& gamma,
+ const size_t& space_order,
+ const size_t& time_order,
+ const CellValue<const double>& rho,
+ const CellValue<const TinyVector<MeshType::Dimension>>& u,
+ const CellValue<const double>& E,
+ const CellValue<const double>& p)
+ : m_dt{dt},
+ m_L{L},
+ m_M{M},
+ m_k{k},
+ m_gamma{gamma},
+ m_spaceOrder{space_order},
+ m_timeOrder{time_order},
+ m_rho{rho},
+ m_u{u},
+ m_E{E},
+ m_p{p},
+ p_mesh{mesh},
+ m_mesh{*mesh}
+ {
+ m_Mj = [&] {
+ CellValue<double> Mj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Mj[cell_id] = m_rho[cell_id] * m_Vj[cell_id]; });
+ return Mj;
+ }();
+
+ m_inv_Mj = [&] {
+ CellValue<double> inv_Mj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Mj[cell_id] = 1. / m_Mj[cell_id]; });
+ return inv_Mj;
+ }();
+
+ m_inv_Vj = [&] {
+ CellValue<double> inv_Vj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Vj[cell_id] = 1. / m_Vj[cell_id]; });
+ return inv_Vj;
+ }();
+ }
+
+ ~EulerKineticSolverLagrangeMultiDLocalOrder2Periodic() = default;
+};
+
+std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiDLocalOrder2Periodic(const double& dt,
+ const size_t& L,
+ const size_t& M,
+ const size_t& k,
+ const double& gamma,
+ const double& eps,
+ const size_t& space_order,
+ const size_t& time_order,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_n,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_n,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_n,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_n)
+{
+ const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({rho_n, u_n, E_n, p_n});
+ if (mesh_v.use_count() == 0) {
+ throw NormalError("rho_n, u_n, E_n and p_n have to be defined on the same mesh");
+ }
+ if (space_order > 1 or time_order > 1) {
+ throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 1");
+ }
+ double eps_bis = eps; // A retirer !!
+ eps_bis += eps_bis; // A retirer !!
+ return std::visit(
+ [&](auto&& p_mesh)
+ -> std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> {
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ if constexpr (MeshType::Dimension != 1) {
+ throw NotImplementedError("periodic bc are only implemented in 1d");
+ } else {
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ auto rho = rho_n->get<DiscreteFunctionP0<const double>>();
+ auto u = u_n->get<DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>>();
+ auto E = E_n->get<DiscreteFunctionP0<const double>>();
+ auto p = p_n->get<DiscreteFunctionP0<const double>>();
+
+ EulerKineticSolverLagrangeMultiDLocalOrder2Periodic solver(p_mesh, dt, L, M, k, gamma, space_order,
+ time_order, rho.cellValues(), u.cellValues(),
+ E.cellValues(), p.cellValues());
+
+ return solver.apply();
+ } else {
+ throw NotImplementedError("Invalid mesh type for multi-D Lagrangian Kinetic solver");
+ }
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.hpp b/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..ef7bb2f0dc4562e9c0738117c17a2d6deb591afe
--- /dev/null
+++ b/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.hpp
@@ -0,0 +1,28 @@
+#ifndef EULER_KINETIC_SOLVER_LAGRANGE_MULTI_D_LOCAL_ORDER_2PERIODIC_HPP
+#define EULER_KINETIC_SOLVER_LAGRANGE_MULTI_D_LOCAL_ORDER_2PERIODIC_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+class IBoundaryConditionDescriptor;
+
+std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiDLocalOrder2Periodic(
+ const double& dt,
+ const size_t& L,
+ const size_t& M,
+ const size_t& k,
+ 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,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+#endif // EULER_KINETIC_SOLVER_LAGRANGE_MULTI_D_LOCAL_ORDER_2PERIODIC_HPP