diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 8e917c929d5216c31411775985720637efe4d526..b82b21a07684be2f4a3bc7470019c547cae8950e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -630,44 +630,41 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("implicit_2_states_euler_o2_rhombus_modified",
- std::function(
+ this->_addBuiltinFunction(
+ "implicit_2_states_euler_o2_rhombus_modified",
+ std::function(
- [](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>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
- const std::shared_ptr<const DiscreteFunctionVariant>& pi,
- const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
- const std::shared_ptr<const DiscreteFunctionVariant>& Cv,
- const std::shared_ptr<const DiscreteFunctionVariant>& entropy,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list,
- const double& dt)
- -> std::tuple<std::shared_ptr<const MeshVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>, double> {
- return ImplicitAcousticO2SolverHandlerModified{ImplicitAcousticO2SolverHandlerModified::
- SolverType::Glace2States,
- rho,
- c,
- u,
- p,
- pi,
- gamma,
- Cv,
- entropy,
- bc_descriptor_list,
- std::vector<std::shared_ptr<
- const IZoneDescriptor>>{},
- dt}
- .solver()
- .apply(dt, rho, u, E, c, p, pi, gamma, Cv, entropy);
- }
+ [](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>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& pi,
+ const std::shared_ptr<const DiscreteFunctionVariant>& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& Cv,
+ const std::shared_ptr<const DiscreteFunctionVariant>& entropy,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, const double& dt)
+ -> std::tuple<std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ double> {
+ return ImplicitAcousticO2SolverHandlerRhombusModified{ImplicitAcousticO2SolverHandlerRhombusModified::
+ SolverType::Glace2States,
+ rho,
+ c,
+ u,
+ p,
+ pi,
+ gamma,
+ Cv,
+ entropy,
+ bc_descriptor_list,
+ std::vector<std::shared_ptr<const IZoneDescriptor>>{},
+ dt}
+ .solver()
+ .apply(dt, rho, u, E, c, p, pi, gamma, Cv, entropy);
+ }
- ));
+ ));
this->_addBuiltinFunction(
"implicit_2_states_euler_o2",
diff --git a/src/scheme/ImplicitAcousticO2SolverRhombus.cpp b/src/scheme/ImplicitAcousticO2SolverRhombus.cpp
index 99e4e867f52fdd4d54f70350ff7d5c7a075cd724..21669540e3c6b654c17427d1ac315c8999b383fa 100644
--- a/src/scheme/ImplicitAcousticO2SolverRhombus.cpp
+++ b/src/scheme/ImplicitAcousticO2SolverRhombus.cpp
@@ -1769,9 +1769,6 @@ class ImplicitAcousticO2SolverHandlerRhombus::ImplicitAcousticO2Solver final
double m = min(new_Vj);
if (m < 0) {
throw NormalError("Negative volume");
- parallel_for(
- mesh->numberOfNodes(),
- PUGS_LAMBDA(const NodeId node_id) { new_xr[node_id] = 0.5 * (new_xr[node_id] + mesh->xr()[node_id]); });
is_positive = false;
}
} while (not is_positive);
diff --git a/src/scheme/ImplicitAcousticO2SolverRhombusModified.cpp b/src/scheme/ImplicitAcousticO2SolverRhombusModified.cpp
index a3df8d3f22680d0da7420b3af513a7e3d8b8c7df..11aa27270a12bca2e6642068d84b10ad7c332da4 100644
--- a/src/scheme/ImplicitAcousticO2SolverRhombusModified.cpp
+++ b/src/scheme/ImplicitAcousticO2SolverRhombusModified.cpp
@@ -1702,243 +1702,250 @@ class ImplicitAcousticO2SolverHandlerRhombusModified::ImplicitAcousticO2Solver f
});
return computed_g_1_exp_S_Cv_inv_g;
}();
- bool CFL = false;
+
double dt_f = dt;
+ double max_tau_error;
+ int number_iter = 0;
do {
- bool test = true;
- double max_tau_error;
- int number_iter = 0;
- do {
- number_iter++;
+ number_iter++;
- NodeValue<const Rd> ur_n = this->_getUr();
- NodeValuePerCell<const Rd> Fjr_n = this->_getFjr(ur_n);
+ NodeValue<const Rd> ur_n = this->_getUr();
+ NodeValuePerCell<const Rd> Fjr_n = this->_getFjr(ur_n);
- this->_computeGradJU_AU(u, p, pi, dt);
+ const auto& cell_to_node_matrix_n = m_mesh.connectivity().cellToNodeMatrix();
- NodeValue<const Rd> ur = this->_getUr();
- NodeValuePerCell<const Rd> Fjr = this->_getFjr(ur);
- // time n+1
+ for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+ double new_tau_fluxes_n = 0;
+ const auto& cell_nodes = cell_to_node_matrix_n[j];
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ new_tau_fluxes_n += dot(m_Djr(j, R), ur_n[r]);
+ }
+ if (-new_tau_fluxes_n > 0) {
+ dt_f = std::min(m_tau[j] * m_Mj[j] / -new_tau_fluxes_n, dt_f);
+ std ::cout << "la c'est " << j << '\n';
+ }
+ }
+ std::cout << "dt = " << dt << "dt_f = " << dt_f << '\n';
+ this->_computeGradJU_AU(u, p, pi, dt_f);
- const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ NodeValue<const Rd> ur = this->_getUr();
+ NodeValuePerCell<const Rd> Fjr = this->_getFjr(ur);
+ // time n+1
- NodeValue<Rd> new_xr = copy(mesh->xr());
- parallel_for(
- mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += 0.5 * dt * (ur[r] + ur_n[r]); });
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
- new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+ NodeValue<Rd> new_xr = copy(mesh->xr());
+ parallel_for(
+ mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += 0.5 * dt_f * (ur[r] + ur_n[r]); });
- CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+ new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
- new_rho = copy(rho.cellValues());
- new_u = copy(u.cellValues());
- new_E = copy(E.cellValues());
+ CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
+ new_rho = copy(rho.cellValues());
+ new_u = copy(u.cellValues());
+ new_E = copy(E.cellValues());
- Rd momentum_fluxes = zero;
- double energy_fluxes = 0;
- for (size_t R = 0; R < cell_nodes.size(); ++R) {
- const NodeId r = cell_nodes[R];
- momentum_fluxes += Fjr(j, R) + Fjr_n(j, R);
- energy_fluxes += dot(Fjr(j, R), ur[r]) + dot(Fjr_n(j, R), ur_n[r]);
- }
- const double dt_over_Mj = dt / m_Mj[j];
- new_u[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
- new_E[j] -= 0.5 * dt_over_Mj * energy_fluxes;
- });
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
- // update Djr
- const NodeValuePerCell<const Rd> new_Cjr = MeshDataManager::instance().getMeshData(*new_mesh).Cjr();
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ momentum_fluxes += Fjr(j, R) + Fjr_n(j, R);
+ energy_fluxes += dot(Fjr(j, R), ur[r]) + dot(Fjr_n(j, R), ur_n[r]);
+ }
+ const double dt_over_Mj = dt_f / m_Mj[j];
+ new_u[j] -= 0.5 * dt_over_Mj * momentum_fluxes;
+ new_E[j] -= 0.5 * dt_over_Mj * energy_fluxes;
+ });
- CellValue<double> new_Vj{m_mesh.connectivity()};
+ // update Djr
+ const NodeValuePerCell<const Rd> new_Cjr = MeshDataManager::instance().getMeshData(*new_mesh).Cjr();
- bool is_positive = true;
- do {
- is_positive = true;
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- new_Vj[j] = 0;
- auto cell_nodes = cell_to_node_matrix[j];
+ CellValue<double> new_Vj{m_mesh.connectivity()};
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- const NodeId node_id = cell_nodes[i_node];
- new_Vj[j] += dot(new_Cjr[j][i_node], new_xr[node_id]);
- }
- new_Vj[j] *= 1. / Dimension;
- });
+ bool is_positive = true;
+ do {
+ is_positive = true;
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ new_Vj[j] = 0;
+ auto cell_nodes = cell_to_node_matrix[j];
- double m = min(new_Vj);
- if (m < 0) {
- // std::cout << "negative volume" << '/n';
- // throw NormalError("Negative volume");
- }
- } while (not is_positive);
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId node_id = cell_nodes[i_node];
+ new_Vj[j] += dot(new_Cjr[j][i_node], new_xr[node_id]);
+ }
+ new_Vj[j] *= 1. / Dimension;
+ });
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+ double m = min(new_Vj);
+ if (m < 0) {
+ std::cout << "la plus petite maille est" << m << '\n';
+ throw NormalError("Negative volume");
+ }
+ } while (not is_positive);
- CellValue<double> new_tau(m_mesh.connectivity());
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
- for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
- double new_tau_fluxes = 0;
- const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t R = 0; R < cell_nodes.size(); ++R) {
- const NodeId r = cell_nodes[R];
- new_tau_fluxes += dot(m_Djr(j, R), ur[r]) + dot(m_Djr(j, R), ur_n[r]);
- }
- new_tau[j] = m_tau[j] + 0.5 * (dt / m_Mj[j]) * new_tau_fluxes;
- // std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
+ CellValue<double> new_tau(m_mesh.connectivity());
+
+ for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+ double new_tau_fluxes = 0;
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R = 0; R < cell_nodes.size(); ++R) {
+ const NodeId r = cell_nodes[R];
+ new_tau_fluxes += dot(m_Djr(j, R), ur[r]) + dot(m_Djr(j, R), ur_n[r]);
}
+ new_tau[j] = m_tau[j] + 0.5 * (dt_f / m_Mj[j]) * new_tau_fluxes;
+ // std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
+ }
- // double sum_tau_rho = 0;
- max_tau_error = 0;
+ // double sum_tau_rho = 0;
+ max_tau_error = 0;
- for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
- // const auto& cell_nodes = cell_to_node_matrix[j];
- // V_j^n+1/tau_j^n+1 - M_j
- // sum_tau_rho += std::abs(new_tau[j] - 1. / new_rho[j]);
- if (m_is_implicit_cell[j]) {
- max_tau_error = std::max(max_tau_error, std::abs(1 - new_tau[j] * new_rho[j]));
- }
+ for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+ // const auto& cell_nodes = cell_to_node_matrix[j];
+ // V_j^n+1/tau_j^n+1 - M_j
+ // sum_tau_rho += std::abs(new_tau[j] - 1. / new_rho[j]);
+ if (m_is_implicit_cell[j]) {
+ max_tau_error = std::max(max_tau_error, std::abs(1 - new_tau[j] * new_rho[j]));
}
- // std::cout << "sum_tau_rho =" << sum_tau_rho << '\n';
- // std::cout << "max_tau_error=" << max_tau_error << '\n';
+ }
+ // std::cout << "sum_tau_rho =" << sum_tau_rho << '\n';
+ // std::cout << "max_tau_error=" << max_tau_error << '\n';
- if constexpr (Dimension == 2) {
- NodeValuePerCell<Rd> Djr{m_mesh.connectivity()};
+ if constexpr (Dimension == 2) {
+ NodeValuePerCell<Rd> Djr{m_mesh.connectivity()};
- parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- const auto& cell_nodes = cell_to_node_matrix[cell_id];
- if (m_is_implicit_cell[cell_id]) {
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- Djr(cell_id, i_node) = 0.5 * (Cjr_n(cell_id, i_node) + new_Cjr(cell_id, i_node));
- }
- } else {
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- Djr(cell_id, i_node) = Cjr_n(cell_id, i_node);
- }
+ parallel_for(
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ if (m_is_implicit_cell[cell_id]) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ Djr(cell_id, i_node) = 0.5 * (Cjr_n(cell_id, i_node) + new_Cjr(cell_id, i_node));
+ }
+ } else {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ Djr(cell_id, i_node) = Cjr_n(cell_id, i_node);
}
- });
+ }
+ });
- m_Djr = Djr;
- } else if constexpr (Dimension == 3) {
- NodeValuePerFace<Rd> Nlr_delta(m_mesh.connectivity());
- const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ m_Djr = Djr;
+ } else if constexpr (Dimension == 3) {
+ NodeValuePerFace<Rd> Nlr_delta(m_mesh.connectivity());
+ const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
- const auto& xr_n = m_mesh.xr();
- const auto& xr_np1 = new_xr;
+ const auto& xr_n = m_mesh.xr();
+ const auto& xr_np1 = new_xr;
- parallel_for(
- m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId l) {
- const auto& face_nodes = face_to_node_matrix[l];
- const size_t nb_nodes = face_nodes.size();
- std::vector<Rd> dxr_n(nb_nodes);
- std::vector<Rd> dxr_np1(nb_nodes);
- for (size_t r = 0; r < nb_nodes; ++r) {
- dxr_n[r] = xr_n[face_nodes[(r + 1) % nb_nodes]] - xr_n[face_nodes[(r + nb_nodes - 1) % nb_nodes]];
- dxr_np1[r] = xr_np1[face_nodes[(r + 1) % nb_nodes]] - xr_np1[face_nodes[(r + nb_nodes - 1) % nb_nodes]];
- }
- const double inv_12_nb_nodes = 1. / (12. * nb_nodes);
- for (size_t r = 0; r < nb_nodes; ++r) {
- Rd Nr = zero;
- for (size_t s = 0; s < nb_nodes; ++s) {
- Nr -= crossProduct((1. / 6) * (2 * (dxr_np1[r] - dxr_n[r]) - (dxr_np1[s] - dxr_n[s])),
- xr_np1[face_nodes[s]] - xr_n[face_nodes[s]]);
- }
- Nr *= inv_12_nb_nodes;
- Nlr_delta(l, r) = Nr;
+ parallel_for(
+ m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId l) {
+ const auto& face_nodes = face_to_node_matrix[l];
+ const size_t nb_nodes = face_nodes.size();
+ std::vector<Rd> dxr_n(nb_nodes);
+ std::vector<Rd> dxr_np1(nb_nodes);
+ for (size_t r = 0; r < nb_nodes; ++r) {
+ dxr_n[r] = xr_n[face_nodes[(r + 1) % nb_nodes]] - xr_n[face_nodes[(r + nb_nodes - 1) % nb_nodes]];
+ dxr_np1[r] = xr_np1[face_nodes[(r + 1) % nb_nodes]] - xr_np1[face_nodes[(r + nb_nodes - 1) % nb_nodes]];
+ }
+ const double inv_12_nb_nodes = 1. / (12. * nb_nodes);
+ for (size_t r = 0; r < nb_nodes; ++r) {
+ Rd Nr = zero;
+ for (size_t s = 0; s < nb_nodes; ++s) {
+ Nr -= crossProduct((1. / 6) * (2 * (dxr_np1[r] - dxr_n[r]) - (dxr_np1[s] - dxr_n[s])),
+ xr_np1[face_nodes[s]] - xr_n[face_nodes[s]]);
}
- });
+ Nr *= inv_12_nb_nodes;
+ Nlr_delta(l, r) = Nr;
+ }
+ });
- const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
- const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+ const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
- NodeValuePerCell<Rd> Djr{m_mesh.connectivity()};
- Djr.fill(zero);
+ NodeValuePerCell<Rd> Djr{m_mesh.connectivity()};
+ Djr.fill(zero);
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
- const auto& cell_faces = cell_to_face_matrix[j];
- const auto& face_is_reversed = cell_face_is_reversed.itemArray(j);
+ const auto& cell_faces = cell_to_face_matrix[j];
+ const auto& face_is_reversed = cell_face_is_reversed.itemArray(j);
- for (size_t L = 0; L < cell_faces.size(); ++L) {
- const FaceId& l = cell_faces[L];
- const auto& face_nodes = face_to_node_matrix[l];
+ for (size_t L = 0; L < cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_node_matrix[l];
- auto local_node_number_in_cell = [&](NodeId node_number) {
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- if (node_number == cell_nodes[i_node]) {
- return i_node;
- }
+ auto local_node_number_in_cell = [&](NodeId node_number) {
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
}
- return std::numeric_limits<size_t>::max();
- };
+ }
+ return std::numeric_limits<size_t>::max();
+ };
- if (face_is_reversed[L]) {
- for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes[rl]);
- Djr(j, R) -= Nlr_delta(l, rl);
- }
- } else {
- for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes[rl]);
- Djr(j, R) += Nlr_delta(l, rl);
- }
+ if (face_is_reversed[L]) {
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Djr(j, R) -= Nlr_delta(l, rl);
+ }
+ } else {
+ for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Djr(j, R) += Nlr_delta(l, rl);
}
}
- });
+ }
+ });
- parallel_for(
- Djr.numberOfValues(), PUGS_LAMBDA(const size_t i) { Djr[i] += 0.5 * (Cjr_n[i] + new_Cjr[i]); });
+ parallel_for(
+ Djr.numberOfValues(), PUGS_LAMBDA(const size_t i) { Djr[i] += 0.5 * (Cjr_n[i] + new_Cjr[i]); });
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- if (not m_is_implicit_cell[cell_id]) {
- const auto& cell_nodes = cell_to_node_matrix[cell_id];
- for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
- Djr(cell_id, i_node) = Cjr_n(cell_id, i_node);
- }
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ if (not m_is_implicit_cell[cell_id]) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ Djr(cell_id, i_node) = Cjr_n(cell_id, i_node);
}
- });
+ }
+ });
- m_Djr = Djr;
+ m_Djr = Djr;
- double max_err = 0;
+ double max_err = 0;
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- if (m_is_implicit_cell[cell_id]) {
- double delta_V = new_Vj[cell_id] - Vj[cell_id];
- const auto& node_list = cell_to_node_matrix[cell_id];
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId node_id = node_list[i_node];
- delta_V -= dt * dot(ur[node_id], Djr[cell_id][i_node]);
- }
- max_err = std::max(max_err, std::abs(delta_V));
+ for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ if (m_is_implicit_cell[cell_id]) {
+ double delta_V = new_Vj[cell_id] - Vj[cell_id];
+ const auto& node_list = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ delta_V -= dt_f * dot(ur[node_id], Djr[cell_id][i_node]);
}
+ max_err = std::max(max_err, std::abs(delta_V));
}
-
- std::cout << rang::fgB::yellow << "Max volume err = " << max_err << rang::fg::reset << '\n';
}
- if constexpr (Dimension > 1) {
- std::cout << rang::fgB::magenta << "number_iter_Djr=" << number_iter << " max_tau_error=" << max_tau_error
- << rang::fg::reset << '\n';
- }
+ std::cout << rang::fgB::yellow << "Max volume err = " << max_err << rang::fg::reset << '\n';
+ }
- // std::cout << "new rho=" << new_rho << '\n';
- } while ((Dimension > 1) and (max_tau_error > 1e-4) and (number_iter < 100));
- if (test) {
- CFL = true;
- } else {
- dt_f = 0.5 * dt_f;
+ if constexpr (Dimension > 1) {
+ std::cout << rang::fgB::magenta << "number_iter_Djr=" << number_iter << " max_tau_error=" << max_tau_error
+ << rang::fg::reset << '\n';
}
- } while (CFL == false);
+
+ // std::cout << "new rho=" << new_rho << '\n';
+ } while ((Dimension > 1) and (max_tau_error > 1e-4) and (number_iter < 100));
// std::cout << "prochaine boucle" << '\n';
implicit_t.pause();
// std::cout << "getA=" << get_A_t << "(" << count_getA << ")"
diff --git a/src/scheme/ImplicitAcousticO2SolverRhombusModified.hpp b/src/scheme/ImplicitAcousticO2SolverRhombusModified.hpp
index 2edf4f8b78546fef972a64b8f795d6cade4c28a1..818a33d383ceb989cef5dc7567ede8a70b37b83f 100644
--- a/src/scheme/ImplicitAcousticO2SolverRhombusModified.hpp
+++ b/src/scheme/ImplicitAcousticO2SolverRhombusModified.hpp
@@ -1,5 +1,5 @@
-#ifndef IMPLICIT_ACOUSTIC_O2_SOLVER_MODIFIED_HPP
-#define IMPLICIT_ACOUSTIC_O2_SOLVER_MODIFIED_HPP
+#ifndef IMPLICIT_ACOUSTIC_O2_SOLVER_RHOMBUS_MODIFIED_HPP
+#define IMPLICIT_ACOUSTIC_O2_SOLVER_RHOMBUS_MODIFIED_HPP
#include <mesh/MeshTraits.hpp>