diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index e6827ea1a20c59946c772509db1a580a9c9e4c1d..437cb5c4ae8d0079b774018de43fe46f924b0ea2 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -2349,144 +2349,6 @@ class EulerKineticSolverMultiDOrder2
});
}
- void
- velocity_limiter(const CellValue<const TinyVector<Dimension>>& u,
- DiscreteFunctionDPk<Dimension, const double>& DPk_p,
- DiscreteFunctionDPk<Dimension, TinyVector<Dimension>>& DPk_uh) const
- {
- StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
- StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
- auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor,
- symmetry_boundary_descriptor_list);
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- auto coefficients_p = DPk_p.coefficients(cell_id);
- TinyVector<Dimension> n;
- for (size_t dim = 1; dim < coefficients_p.size(); ++dim) {
- n[dim - 1] = coefficients_p[dim];
- }
- if (dot(n, n) != 0) {
- n = 1. / sqrt(dot(n, n)) * n;
-
- TinyVector<Dimension> v;
- bool is_colinear = false;
- v[0] = 1;
- for (size_t dim = 1; dim < Dimension; ++dim) {
- v[dim] = 0;
- if (n[0] != 0 and n[dim] == 0) {
- is_colinear = true;
- } else {
- is_colinear = false;
- }
- }
- if (is_colinear) {
- for (size_t dim = 1; dim < Dimension; ++dim) {
- if (dim != 1) {
- v[dim] = 0;
- } else {
- v[dim] = 1;
- }
- }
- }
-
- TinyVector<Dimension> t = v - dot(v, n) * n;
- t = 1. / sqrt(dot(t, t)) * t;
- TinyVector<Dimension> l = zero;
-
- if constexpr (Dimension == 3) {
- l[0] = n[1] * t[2] - n[2] * t[1];
- l[1] = n[2] * t[0] - n[0] * t[2];
- l[0] = n[0] * t[1] - n[1] * t[0];
- l = 1. / sqrt(dot(l, l)) * l;
- }
-
- const double un = dot(u[cell_id], n);
- const double ut = dot(u[cell_id], t);
- const double ul = dot(u[cell_id], l);
-
- double un_min = un;
- double un_max = un;
- double ut_min = ut;
- double ut_max = ut;
- double ul_min = ul;
- double ul_max = ul;
-
- const auto cell_stencil = stencil[cell_id];
- for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
- un_min = std::min(un_min, dot(u[cell_stencil[i_cell]], n));
- un_max = std::max(un_max, dot(u[cell_stencil[i_cell]], n));
- ut_min = std::min(ut_min, dot(u[cell_stencil[i_cell]], t));
- ut_max = std::max(ut_max, dot(u[cell_stencil[i_cell]], t));
- ul_min = std::min(ul_min, dot(u[cell_stencil[i_cell]], l));
- ul_max = std::max(ul_max, dot(u[cell_stencil[i_cell]], l));
- }
-
- const double eps = 1E-14;
- double coef_n = 1;
- double lambda_n = 1.;
- double coef_t = 1;
- double lambda_t = 1.;
- double coef_l = 1;
- double lambda_l = 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 un_bar_i = dot(DPk_uh[cell_id](m_xj[cell_i_id]), n);
- const double Delta_n = (un_bar_i - un);
- const double ut_bar_i = dot(DPk_uh[cell_id](m_xj[cell_i_id]), t);
- const double Delta_t = (ut_bar_i - ut);
- const double ul_bar_i = dot(DPk_uh[cell_id](m_xj[cell_i_id]), l);
- const double Delta_l = (ul_bar_i - ul);
-
- if (Delta_n > eps) {
- coef_n = std::min(1., (un_max - un) / Delta_n);
- } else if (Delta_n < -eps) {
- coef_n = std::min(1., (un_min - un) / Delta_n);
- }
- lambda_n = std::min(lambda_n, coef_n);
-
- if (Delta_t > eps) {
- coef_t = std::min(1., (ut_max - ut) / Delta_t);
- } else if (Delta_t < -eps) {
- coef_t = std::min(1., (ut_min - ut) / Delta_t);
- }
- lambda_t = std::min(lambda_t, coef_t);
-
- if (Delta_l > eps) {
- coef_l = std::min(1., (ul_max - ul) / Delta_l);
- } else if (Delta_l < -eps) {
- coef_l = std::min(1., (ul_min - ul) / Delta_l);
- }
- lambda_l = std::min(lambda_l, coef_l);
- }
-
- auto coefficients = DPk_uh.coefficients(cell_id);
- coefficients[0] = (1 - lambda_n) * dot(u[cell_id], n) * n + lambda_n * dot(coefficients[0], n) * n +
- (1 - lambda_t) * dot(u[cell_id], t) * t + lambda_t * dot(coefficients[0], t) * t +
- (1 - lambda_l) * dot(u[cell_id], l) * l + lambda_l * dot(coefficients[0], l) * l;
- for (size_t i = 1; i < coefficients.size(); ++i) {
- coefficients[i] = lambda_n * dot(coefficients[i], n) * n + lambda_t * dot(coefficients[i], t) * t +
- lambda_l * dot(coefficients[i], l) * l;
- }
- }
- });
- }
-
- const CellValue<const double>
- build_pressure(const CellValue<const double> rho,
- const CellValue<const TinyVector<Dimension>> rho_u,
- const CellValue<const double> rho_E)
- {
- CellValue<double> p{p_mesh->connectivity()};
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
- p[cell_id] = (m_gamma - 1) * rho_e;
- });
- return p;
- }
-
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -2647,15 +2509,13 @@ class EulerKineticSolverMultiDOrder2
DiscreteFunctionDPk<Dimension, double> rho_E_lim = copy(DPk_rho_E);
if (limitation) {
- const CellValue<const double> p = build_pressure(rho, rho_u, rho_E);
- auto reconstruction_pressure =
- PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, p));
- auto DPk_p = reconstruction_pressure[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
-
scalar_limiter_equilibrium(rho, rho_lim);
scalar_limiter_equilibrium(rho_E, rho_E_lim);
vector_limiter_equilibrium(rho_u, rho_u_lim);
- // velocity_limiter(rho_u, DPk_p, rho_u_lim);
+
+ synchronize(rho_lim.cellArrays());
+ synchronize(rho_u_lim.cellArrays());
+ synchronize(rho_E_lim.cellArrays());
}
if (m_scheme == 1) {
@@ -2752,9 +2612,7 @@ class EulerKineticSolverMultiDOrder2
NodeArray<double> Fr_rho = compute_Flux_glace<double>(Fn_rho_lim, is_boundary_node);
NodeArray<TinyVector<Dimension>> Fr_rho_u =
compute_Flux_glace<TinyVector<Dimension>>(Fn_rho_u_lim, is_boundary_node);
- NodeArray<double>
-
- Fr_rho_E = compute_Flux_glace<double>(Fn_rho_E_lim, is_boundary_node);
+ NodeArray<double> Fr_rho_E = compute_Flux_glace<double>(Fn_rho_E_lim, is_boundary_node);
apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list);
@@ -2803,15 +2661,13 @@ class EulerKineticSolverMultiDOrder2
DiscreteFunctionDPk<Dimension, double> rho_E_np1_lim = copy(DPk_rho_E_np1);
if (limitation) {
- const CellValue<const double> p_np1 = build_pressure(rho_np1, rho_u_np1, rho_E_np1);
- auto reconstruction_pressure_np1 =
- PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0(p_mesh, p_np1));
- auto DPk_p_np1 = reconstruction_pressure_np1[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
-
scalar_limiter_equilibrium(rho_np1, rho_np1_lim);
scalar_limiter_equilibrium(rho_E_np1, rho_E_np1_lim);
vector_limiter_equilibrium(rho_u_np1, rho_u_np1_lim);
- // velocity_limiter(rho_u_np1, DPk_p_np1, rho_u_np1_lim);
+
+ synchronize(rho_np1_lim.cellArrays());
+ synchronize(rho_u_np1_lim.cellArrays());
+ synchronize(rho_E_np1_lim.cellArrays());
}
if (m_scheme == 1) {
@@ -2840,6 +2696,7 @@ class EulerKineticSolverMultiDOrder2
deltaLFnp1_rho = compute_deltaLFn_face<double>(Fl_rho_np1);
deltaLFnp1_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1);
deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1);
+
} else if (m_scheme == 2) {
NodeArray<double> Fr_rho_np1{p_mesh->connectivity(), nb_waves};
NodeArray<TinyVector<Dimension>> Fr_rho_u_np1{p_mesh->connectivity(), nb_waves};