diff --git a/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp b/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp
index 6c188eec384938f2a80d377a9c5c9d04da21e0d5..5a3164385581a51460f0213d3a4c41e7d09a6c2b 100644
--- a/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp
+++ b/src/scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.cpp
@@ -771,61 +771,96 @@ class EulerKineticSolverAcoustic2LagrangeFVOrder2
NodeArray<double> p_flux(m_mesh.connectivity(), m_timeOrder);
CellValue<double> new_p{m_mesh.connectivity()};
+ CellValue<TinyVector<Dimension>> new_u{m_mesh.connectivity()};
+ std::shared_ptr<const MeshType> new_mesh;
- for (size_t i_time = 0; i_time < m_timeOrder; ++i_time) {
- if (i_time == 0) {
- vec_F = compute_M(m_p, m_u);
- } else {
- vec_F = compute_M(p, u);
- }
- F_p = get<0>(vec_F);
- F_u = get<1>(vec_F);
- F0_p = copy(F_p);
- F0_u = copy(F_u);
-
- NodeArray<double> Flux_Fn_p(m_mesh.connectivity(), nb_waves);
- NodeArray<double> Flux_Fn_u(m_mesh.connectivity(), nb_waves);
-
- if (m_spaceOrder == 1) {
- Flux_Fn_p = compute_Flux1(F_p);
- Flux_Fn_u = compute_Flux1(F_u);
- } else if (m_spaceOrder == 2) {
- Flux_Fn_p = compute_Polynomial_reconstruction_2(F_p);
- Flux_Fn_u = compute_Polynomial_reconstruction_2(F_u);
- } else if (m_spaceOrder == 3) {
- Flux_Fn_p = compute_Polynomial_reconstruction_3(F_p);
- Flux_Fn_u = compute_Polynomial_reconstruction_3(F_u);
- } else {
- throw NotImplementedError("Only implemented for space orders 1, 2 and 3");
- }
+ while (Mood_activation) {
+ Mood_activation = false;
- if (m_TVDLimitation == 1) {
- Flux_Fn_p = TVD(Flux_Fn_p, F0_p, c[i_time]);
- Flux_Fn_u = TVD(Flux_Fn_u, F0_u, c[i_time]);
- }
+ for (size_t i_time = 0; i_time < m_timeOrder; ++i_time) {
+ if (i_time == 0) {
+ vec_F = compute_M(m_p, m_u);
+ } else {
+ vec_F = compute_M(p, u);
+ }
+ F_p = get<0>(vec_F);
+ F_u = get<1>(vec_F);
+ F0_p = copy(F_p);
+ F0_u = copy(F_u);
+
+ NodeArray<double> Flux_Fn_p(m_mesh.connectivity(), nb_waves);
+ NodeArray<double> Flux_Fn_u(m_mesh.connectivity(), nb_waves);
+ NodeArray<double> Flux1_Fn_p = compute_Flux1(F_p);
+ NodeArray<double> Flux1_Fn_u = compute_Flux1(F_u);
+
+ if (m_spaceOrder == 1) {
+ Flux_Fn_p = compute_Flux1(F_p);
+ Flux_Fn_u = compute_Flux1(F_u);
+ } else if (m_spaceOrder == 2) {
+ Flux_Fn_p = compute_Polynomial_reconstruction_2(F_p);
+ Flux_Fn_u = compute_Polynomial_reconstruction_2(F_u);
+ } else if (m_spaceOrder == 3) {
+ Flux_Fn_p = compute_Polynomial_reconstruction_3(F_p);
+ Flux_Fn_u = compute_Polynomial_reconstruction_3(F_u);
+ } else {
+ throw NotImplementedError("Only implemented for space orders 1, 2 and 3");
+ }
- if (i_time == 0) {
- for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
+ if (m_MOODLimitation) {
for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
- au_flux[node_id][j_time][0] = 0;
- p_flux[node_id][j_time] = 0;
+ if (not node_is_valid[node_id]) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ Flux_Fn_p[node_id][i_wave] = Flux1_Fn_p[node_id][i_wave];
+ Flux_Fn_u[node_id][i_wave] = Flux1_Fn_u[node_id][i_wave];
+ }
+ }
+ }
+ }
+
+ if (m_TVDLimitation == 1) {
+ Flux_Fn_p = TVD(Flux_Fn_p, F0_p, c[i_time]);
+ Flux_Fn_u = TVD(Flux_Fn_u, F0_u, c[i_time]);
+ }
+
+ if (i_time == 0) {
+ for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
+ for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
+ au_flux[node_id][j_time][0] = 0;
+ p_flux[node_id][j_time] = 0;
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ au_flux[node_id][j_time][0] += m_lambda_vector[i_wave] * Flux_Fn_p[node_id][i_wave];
+ p_flux[node_id][j_time] += m_lambda_vector[i_wave] * Flux_Fn_u[node_id][i_wave];
+ }
+ u_flux[node_id][j_time][0] = au_flux[node_id][j_time][0] / (lambda * lambda);
+ }
+ }
+
+ } else {
+ for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
+ au_flux[node_id][i_time][0] = 0;
+ p_flux[node_id][i_time] = 0;
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
- au_flux[node_id][j_time][0] += m_lambda_vector[i_wave] * Flux_Fn_p[node_id][i_wave];
- p_flux[node_id][j_time] += m_lambda_vector[i_wave] * Flux_Fn_u[node_id][i_wave];
+ au_flux[node_id][i_time][0] += m_lambda_vector[i_wave] * Flux_Fn_p[node_id][i_wave];
+ p_flux[node_id][i_time] += m_lambda_vector[i_wave] * Flux_Fn_u[node_id][i_wave];
}
- u_flux[node_id][j_time][0] = au_flux[node_id][j_time][0] / (lambda * lambda);
+ u_flux[node_id][i_time][0] = au_flux[node_id][i_time][0] / (lambda * lambda);
}
}
- } else {
- for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
- au_flux[node_id][i_time][0] = 0;
- p_flux[node_id][i_time] = 0;
- for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
- au_flux[node_id][i_time][0] += m_lambda_vector[i_wave] * Flux_Fn_p[node_id][i_wave];
- p_flux[node_id][i_time] += m_lambda_vector[i_wave] * Flux_Fn_u[node_id][i_wave];
+ for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ const double dt_over_Mj = (m_dt * m_inv_Mj[cell_id]);
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ const NodeId left_node_id = cell_nodes[0];
+ const NodeId right_node_id = cell_nodes[1];
+ u[cell_id][0] = m_u[cell_id][0];
+ tau[cell_id] = 1. / m_rho[cell_id];
+ for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
+ u[cell_id][0] -= B[j_time] * dt_over_Mj * (p_flux[right_node_id][j_time] - p_flux[left_node_id][j_time]);
+
+ tau[cell_id] +=
+ B[j_time] * dt_over_Mj * (u_flux[right_node_id][j_time][0] - u_flux[left_node_id][j_time][0]);
}
- u_flux[node_id][i_time][0] = au_flux[node_id][i_time][0] / (lambda * lambda);
+ p[cell_id] = m_p[cell_id] / std::pow(m_rho[cell_id] * tau[cell_id], m_gamma);
}
}
@@ -834,55 +869,70 @@ class EulerKineticSolverAcoustic2LagrangeFVOrder2
const auto& cell_nodes = cell_to_node_matrix[cell_id];
const NodeId left_node_id = cell_nodes[0];
const NodeId right_node_id = cell_nodes[1];
- u[cell_id][0] = m_u[cell_id][0];
- tau[cell_id] = 1. / m_rho[cell_id];
+
+ u[cell_id][0] = m_u[cell_id][0];
+ E[cell_id] = m_E[cell_id];
for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
u[cell_id][0] -= B[j_time] * dt_over_Mj * (p_flux[right_node_id][j_time] - p_flux[left_node_id][j_time]);
- tau[cell_id] += B[j_time] * dt_over_Mj * (u_flux[right_node_id][j_time][0] - u_flux[left_node_id][j_time][0]);
+ E[cell_id] -= B[j_time] * dt_over_Mj *
+ (p_flux[right_node_id][j_time] * u_flux[right_node_id][j_time][0] -
+ p_flux[left_node_id][j_time] * u_flux[left_node_id][j_time][0]);
}
- p[cell_id] = m_p[cell_id] / std::pow(m_rho[cell_id] * tau[cell_id], m_gamma);
}
- }
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- const double dt_over_Mj = (m_dt * m_inv_Mj[cell_id]);
- const auto& cell_nodes = cell_to_node_matrix[cell_id];
- const NodeId left_node_id = cell_nodes[0];
- const NodeId right_node_id = cell_nodes[1];
+ NodeValue<TinyVector<Dimension>> new_xr = copy(m_mesh.xr());
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
+ new_xr[node_id][0] += B[j_time] * m_dt * u_flux[node_id][j_time][0];
+ }
+ });
- u[cell_id][0] = m_u[cell_id][0];
- E[cell_id] = m_E[cell_id];
- for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
- u[cell_id][0] -= B[j_time] * dt_over_Mj * (p_flux[right_node_id][j_time] - p_flux[left_node_id][j_time]);
+ new_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), new_xr);
- E[cell_id] -= B[j_time] * dt_over_Mj *
- (p_flux[right_node_id][j_time] * u_flux[right_node_id][j_time][0] -
- p_flux[left_node_id][j_time] * u_flux[left_node_id][j_time][0]);
- }
- }
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
- NodeValue<TinyVector<Dimension>> new_xr = copy(m_mesh.xr());
- parallel_for(
- m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- for (size_t j_time = 0; j_time < m_timeOrder; ++j_time) {
- new_xr[node_id][0] += B[j_time] * m_dt * u_flux[node_id][j_time][0];
- }
- });
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rho[cell_id] = m_Mj[cell_id] / new_Vj[cell_id]; });
- std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), new_xr);
+ new_u = copy(u.cellValues());
- CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ new_p[cell_id] = (m_gamma - 1) * rho[cell_id] * (E[cell_id] - 0.5 * new_u[cell_id][0] * new_u[cell_id][0]);
+ });
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rho[cell_id] = m_Mj[cell_id] / new_Vj[cell_id]; });
+ if (m_MOODLimitation) {
+ for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ const auto& cell_to_node = cell_to_node_matrix[cell_id];
+ NodeId left_node = cell_to_node[0];
+ NodeId right_node = cell_to_node[1];
+ const size_t N = m_mesh.numberOfCells();
- CellValue<TinyVector<Dimension>> new_u = copy(u.cellValues());
+ const CellId prev_cell_id = (cell_id + N - 1) % N;
+ const CellId next_cell_id = (cell_id + 1) % N;
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- new_p[cell_id] = (m_gamma - 1) * rho[cell_id] * (E[cell_id] - 0.5 * new_u[cell_id][0] * new_u[cell_id][0]);
- });
+ if (rho[cell_id] <= 0) {
+ Mood_activation = true;
+ cell_is_valid[cell_id] = false;
+ cell_is_valid[prev_cell_id] = false;
+ cell_is_valid[next_cell_id] = false;
+ node_is_valid[left_node] = false;
+ node_is_valid[right_node] = false;
+ } else {
+ if (new_p[cell_id] <= 0) {
+ Mood_activation = true;
+ cell_is_valid[cell_id] = false;
+ cell_is_valid[prev_cell_id] = false;
+ cell_is_valid[next_cell_id] = false;
+ node_is_valid[left_node] = false;
+ node_is_valid[right_node] = false;
+ }
+ }
+ }
+ }
+ }
return std::make_tuple(std::make_shared<MeshVariant>(new_mesh),
std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<const double>(new_mesh, rho)),
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2.cpp b/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2.cpp
index 260520c197989e9e8639220f6f64f437316c5274..a639a0a342744dd1f2d9331685c67f098e1f011b 100644
--- a/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2.cpp
+++ b/src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2.cpp
@@ -1062,8 +1062,8 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2
if (limitation) {
const CellValue<const TinyMatrix<Dimension>> grad_u = compute_grad_u(u_lim);
scalar_limiter(p, p_lim);
- vector_limiter_symmetry(u, u_lim, grad_u);
- // vector_limiter(u, u_lim);
+ // vector_limiter_symmetry(u, u_lim, grad_u);
+ vector_limiter(u, u_lim);
synchronize(p_lim.cellArrays());
synchronize(u_lim.cellArrays());
@@ -1124,8 +1124,8 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2
if (limitation) {
const CellValue<const TinyMatrix<Dimension>> grad_u = compute_grad_u(u_star_lim);
scalar_limiter(pj_star, p_star_lim);
- vector_limiter_symmetry(uj_star, u_star_lim, grad_u);
- // vector_limiter(uj_star, u_star_lim);
+ // vector_limiter_symmetry(uj_star, u_star_lim, grad_u);
+ vector_limiter(uj_star, u_star_lim);
synchronize(p_star_lim.cellArrays());
synchronize(u_star_lim.cellArrays());