From 062d9e1db8ff31997e96f14e09e5437b303cde53 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 20 Oct 2023 19:03:37 +0200
Subject: [PATCH] Few fixes and factorization
---
src/scheme/ImplicitAcousticSolver.cpp | 494 +++++++++++++++-----------
1 file changed, 282 insertions(+), 212 deletions(-)
diff --git a/src/scheme/ImplicitAcousticSolver.cpp b/src/scheme/ImplicitAcousticSolver.cpp
index 2b74689bf..7cfed8708 100644
--- a/src/scheme/ImplicitAcousticSolver.cpp
+++ b/src/scheme/ImplicitAcousticSolver.cpp
@@ -106,6 +106,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CellValue<const double> m_g_1_exp_S_Cv_inv_g;
CellValue<const double> m_tau_iter;
+ CellValue<const Rd> m_u;
CellValue<const double> m_tau;
CellValue<const double> m_Mj;
@@ -358,15 +359,15 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
auto is_boundary_node = m_mesh.connectivity().isBoundaryNode();
non_zeros.fill(0);
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- if (m_is_implicit_cell[cell_id]) {
- const auto& cell_node = cell_to_node_matrix[cell_id];
+ for (CellId cell_j_id = 0; cell_j_id < m_mesh.numberOfCells(); ++cell_j_id) {
+ if (m_is_implicit_cell[cell_j_id]) {
+ const auto& cell_j_nodes = cell_to_node_matrix[cell_j_id];
std::vector<size_t> nb_neighbors;
nb_neighbors.reserve(30);
- for (size_t id_node = 0; id_node < cell_node.size(); ++id_node) {
- if (not is_boundary_node[cell_node[id_node]]) {
- NodeId node_id = cell_node[id_node];
+ for (size_t id_node = 0; id_node < cell_j_nodes.size(); ++id_node) {
+ if (not is_boundary_node[cell_j_nodes[id_node]]) {
+ NodeId node_id = cell_j_nodes[id_node];
const auto& node_cell = node_to_cell_matrix[node_id];
for (size_t cell_i = 0; cell_i < node_cell.size(); ++cell_i) {
CellId id_cell_i = node_cell[cell_i];
@@ -380,11 +381,11 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
auto last = std::unique(nb_neighbors.begin(), nb_neighbors.end());
nb_neighbors.resize(std::distance(nb_neighbors.begin(), last));
- size_t line_index_p = mapP(cell_id);
+ size_t line_index_p = mapP(cell_j_id);
non_zeros[line_index_p] = Dimension * nb_neighbors.size() + 1;
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- size_t line_index_u = mapU(i_dimension, cell_id);
+ size_t line_index_u = mapU(i_dimension, cell_j_id);
non_zeros[line_index_u] = nb_neighbors.size() + 1;
}
@@ -400,21 +401,22 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const auto& node_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
- for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- for (size_t cell_i = 0; cell_i < node_cell.size(); ++cell_i) {
- const CellId id_cell_i = node_cell[cell_i];
- if (m_is_implicit_cell[id_cell_i]) {
- const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
- const int i_index_u = mapU(i_dimension, id_cell_i);
- for (size_t cell_j = 0; cell_j < node_cell.size(); ++cell_j) {
- const CellId id_cell_j = node_cell[cell_j];
- if (m_is_implicit_cell[id_cell_j]) {
- const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
- const int j_index_p = mapP(id_cell_j);
-
- const Rd coef = m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] * m_Djr(id_cell_j, node_nb_in_j);
- A(j_index_p, i_index_u) += coef[i_dimension];
- A(i_index_u, j_index_p) -= coef[i_dimension];
+ for (size_t cell_i = 0; cell_i < node_cell.size(); ++cell_i) {
+ const CellId id_cell_i = node_cell[cell_i];
+ if (m_is_implicit_cell[id_cell_i]) {
+ const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
+ for (size_t cell_j = 0; cell_j < node_cell.size(); ++cell_j) {
+ const CellId id_cell_j = node_cell[cell_j];
+ if (m_is_implicit_cell[id_cell_j]) {
+ const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
+ const int j_index_p = mapP(id_cell_j);
+
+ const Rd Bji = m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] * m_Djr(id_cell_j, node_nb_in_j);
+ for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ const int i_index_u = mapU(i_dimension, id_cell_i);
+
+ A(j_index_p, i_index_u) += Bji[i_dimension];
+ A(i_index_u, j_index_p) -= Bji[i_dimension];
}
}
}
@@ -584,7 +586,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
if (m_is_implicit_cell[j]) {
size_t k = mapP(j);
- computed_tau_iter[j] = m_g_1_exp_S_Cv_inv_g[j] * pow(-Uk[k] + pi[j], -m_inv_gamma[j]);
+ computed_tau_iter[j] = m_g_1_exp_S_Cv_inv_g[j] * std::pow(-Uk[k] + pi[j], -m_inv_gamma[j]);
}
});
return computed_tau_iter;
@@ -899,6 +901,80 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
return F;
}
+ Vector<double>
+ _getF_new(const CRSMatrix<double, int>& A,
+ const Vector<double>& Un,
+ const Vector<double>& Uk,
+ const DiscreteScalarFunction& p,
+ const DiscreteScalarFunction& pi,
+ const DiscreteVectorFunction& u,
+ const double dt)
+ {
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+ NodeValue<const Rd> ur = this->_getUr();
+ NodeValuePerCell<const Rd> Fjr = this->_getFjr(ur);
+
+ m_tau_iter = [&]() {
+ CellValue<double> computed_tau_iter(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ if (m_is_implicit_cell[j]) {
+ size_t k = mapP(j);
+ computed_tau_iter[j] = m_g_1_exp_S_Cv_inv_g[j] * std::pow(-Uk[k] + pi[j], -m_inv_gamma[j]);
+ }
+ });
+ return computed_tau_iter;
+ }();
+
+ CellValue<double> volume_fluxes_sum{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ double sum = 0;
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ const NodeId node_id = cell_nodes[i_node];
+ sum += dot(m_Djr(cell_id, i_node), ur[node_id]);
+ }
+ volume_fluxes_sum[cell_id] = sum;
+ });
+
+ CellValue<Rd> momentum_fluxes_sum{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ Rd sum = zero;
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+ for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+ sum += Fjr(cell_id, i_node);
+ }
+ momentum_fluxes_sum[cell_id] = sum;
+ });
+
+ Vector<double> F{Un.size()};
+ F.fill(0);
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ if (m_is_implicit_cell[cell_id]) {
+ const size_t pj_index = mapP(cell_id);
+ F[pj_index] =
+ m_Mj[cell_id] / dt *
+ (m_g_1_exp_S_Cv_inv_g[cell_id] * std::pow(-Uk[pj_index] + pi[cell_id], -m_inv_gamma[cell_id]) -
+ m_tau[cell_id]) -
+ volume_fluxes_sum[cell_id];
+
+ for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ const size_t uj_i_index = mapU(i_dimension, cell_id);
+
+ F[uj_i_index] =
+ m_Mj[cell_id] / dt * (Uk[uj_i_index] - Un[uj_i_index]) + momentum_fluxes_sum[cell_id][i_dimension];
+ }
+ }
+ });
+
+ return F;
+ }
+
CRSMatrixDescriptor<double>
_getHessianJ(const double dt, const DiscreteScalarFunction& pi, const Vector<double>& Uk)
{
@@ -948,9 +1024,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
(Dimension + 1) * m_number_of_implicit_cells, non_zeros};
const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
-
- const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
count_getGradF++;
static bool getGradF_is_started = false;
@@ -986,14 +1059,14 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const size_t node_nb_in_j_cell = node_local_number_in_its_cells[j_cell];
const int j_index_p = mapP(cell_id_j);
- const auto invArCir = inv_Ar * Cjr(cell_id_j, node_nb_in_j_cell);
+ const auto invArDjr = inv_Ar * m_Djr(cell_id_j, node_nb_in_j_cell);
for (size_t i_cell = 0; i_cell < node_cell.size(); ++i_cell) {
const CellId cell_id_i = node_cell[i_cell];
if (m_is_implicit_cell[cell_id_i]) {
const size_t node_nb_in_cell_i = node_local_number_in_its_cells[i_cell];
const int i_index_p = mapP(cell_id_i);
- Hess_J(j_index_p, i_index_p) += dot(invArCir, Cjr(cell_id_i, node_nb_in_cell_i));
+ Hess_J(j_index_p, i_index_p) += dot(invArDjr, m_Djr(cell_id_i, node_nb_in_cell_i));
}
}
}
@@ -1142,7 +1215,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
const int j_index_p = mapP(id_cell_j);
Hess_J(i_index_p, j_index_p) +=
- dot(Cjr(id_cell_i, node_nb_in_i), inverse_ArxP[node_id] * Cjr(id_cell_j, node_nb_in_j));
+ dot(m_Djr(id_cell_i, node_nb_in_i), inverse_ArxP[node_id] * m_Djr(id_cell_j, node_nb_in_j));
}
}
}
@@ -1203,7 +1276,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// const int j_index_p = mapP(id_cell_j);
// Hess_J(i_index_p, j_index_p) +=
- // dot(Cjr(id_cell_i, node_nb_in_i), m_inv_Ar[node_id] * Cjr(id_cell_j, node_nb_in_j));
+ // dot(m_Djr(id_cell_i, node_nb_in_i), m_inv_Ar[node_id] * m_Djr(id_cell_j, node_nb_in_j));
// }
// }
// }
@@ -1428,7 +1501,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_number_of_implicit_cells > 0) {
do {
nb_iter++;
- Vector<double> f = this->_getF(A, Un, Uk, p, pi, u, dt);
+ // Vector<double> f_old = this->_getF(A, Un, Uk, p, pi, u, dt);
+ Vector<double> f = this->_getF_new(A, Un, Uk, p, pi, u, dt);
std::cout << "building gradf: ";
CRSMatrix<double> gradient_f = this->_getGradientF(A, Uk, pi, dt);
@@ -1466,7 +1540,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
std::cout << rang::style::bold << "norm resid = " << l2Norm(f - gradient_f * sol) << rang::style::reset << '\n';
// std::exit(0);
- Vector<double> U_next = Uk - sol;
+ Vector<double> U_next = Uk - 0.2 * sol;
Array<const double> abs_sol = [&]() {
Array<double> compute_abs_sol{sol.size()};
@@ -1477,24 +1551,24 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
norm_inf_sol = max(abs_sol);
- Uk = U_next;
-
std::cout << "nb_iter= " << nb_iter << " | ratio Newton = " << norm_inf_sol / norm_inf_Un << "\n";
// when it is a hard case we relax newton
size_t neg_pressure_count = 0;
double min_pressure = 0;
- for (CellId j = 0; j < m_number_of_implicit_cells; ++j) {
- size_t k = mapP(j);
- if (-U_next[k] + pi[j] < 0) {
+ for (CellId cell_id = 0; cell_id < m_number_of_implicit_cells; ++cell_id) {
+ size_t k = mapP(cell_id);
+ if (-U_next[k] + pi[cell_id] < 0) {
+ std::cout << " neg p: cell_id=" << cell_id << '\n';
++neg_pressure_count;
- min_pressure = std::min(min_pressure, -U_next[k] + pi[j]);
- U_next[k] = Uk[k] + sol[k];
+ min_pressure = std::min(min_pressure, -U_next[k] + pi[cell_id]);
+ U_next[k] = Uk[k];
}
}
if (neg_pressure_count > 0) {
std::cout << rang::fgB::magenta << "p est negatif sur " << neg_pressure_count
<< " mailles min=" << min_pressure << rang::fg::reset << '\n';
+ std::exit(0);
}
// Uk -= sol;
// if (neg_pressure) {
@@ -1504,33 +1578,36 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// Uk -= sol;
// }
- } while ((norm_inf_sol > 1e-15 * norm_inf_Un) and (nb_iter < 200));
- }
- m_predicted_u = [&] {
- CellValue<Rd> predicted_u = copy(u.cellValues());
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- if (m_is_implicit_cell[cell_id]) {
- Rd vector_u = zero;
- for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- vector_u[i_dimension] = Uk[mapU(i_dimension, cell_id)];
+ Uk = U_next;
+
+ m_predicted_u = [&] {
+ CellValue<Rd> predicted_u = copy(u.cellValues());
+ for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ if (m_is_implicit_cell[cell_id]) {
+ Rd vector_u = zero;
+ for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ vector_u[i_dimension] = Uk[mapU(i_dimension, cell_id)];
+ }
+ predicted_u[cell_id] = vector_u;
+ }
+ // std::cout << "u[" << cell_id << "]=" << predicted_u[cell_id] << '\n';
}
- predicted_u[cell_id] = vector_u;
- }
- // std::cout << "u[" << cell_id << "]=" << predicted_u[cell_id] << '\n';
- }
- return predicted_u;
- }();
+ return predicted_u;
+ }();
- m_predicted_p = [&] {
- CellValue<double> predicted_p = copy(p.cellValues());
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- if (m_is_implicit_cell[cell_id]) {
- predicted_p[cell_id] = -Uk[mapP(cell_id)];
- }
- // std::cout << "p[" << cell_id << "]=" << predicted_p[cell_id] << '\n';
- }
- return predicted_p;
- }();
+ m_predicted_p = [&] {
+ CellValue<double> predicted_p = copy(p.cellValues());
+ for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ if (m_is_implicit_cell[cell_id]) {
+ predicted_p[cell_id] = -Uk[mapP(cell_id)];
+ }
+ // std::cout << "p[" << cell_id << "]=" << predicted_p[cell_id] << '\n';
+ }
+ return predicted_p;
+ }();
+
+ } while ((norm_inf_sol > 1e-15 * norm_inf_Un) or (nb_iter < 100));
+ }
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
// faudrait mettre p+pi
@@ -1538,11 +1615,138 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
std::cout << "pression negative pour la maille" << j << '\n';
}
}
+ }
+
+ NodeValue<const Rd>
+ _getUr() const
+ {
+ const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ NodeValue<Rd> b{m_mesh.connectivity()};
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ Rd br = zero;
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ br += m_Ajr(J, R) * m_predicted_u[J] + m_predicted_p[J] * m_Djr(J, R);
+ }
+
+ b[r] = br;
+ });
+
+ for (const auto& boundary_condition : m_boundary_condition_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
+ if constexpr (Dimension == 1) {
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+ const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ const auto& node_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ Assert(node_cell_list.size() == 1);
+
+ CellId node_cell_id = node_cell_list[0];
+ size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+ b[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+ });
+ } else {
+ const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+ const auto& face_to_cell_matrix = m_mesh.connectivity().faceToCellMatrix();
+ const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ const auto& face_local_numbers_in_their_cells = m_mesh.connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_cell_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& value_list = bc.valueList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+ NodeId node_id = face_nodes[i_node];
+ b[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ // const auto cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ auto is_boundary_node = m_mesh.connectivity().isBoundaryNode();
+
+ const Rd& n = bc.outgoingNormal();
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+ const Array<const NodeId>& node_list = bc.nodeList();
+ parallel_for(
+ bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+ const NodeId r = node_list[r_number];
+ b[r] = P * b[r];
+ });
+
+ } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ const auto& value = value_list[i_node];
+ b[node_id] = value;
+ });
+ } else {
+ throw UnexpectedError("boundary condition not handled 5");
+ }
+ },
+ boundary_condition);
+ }
+
+ const NodeValue<Rd> computed_ur(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = m_inv_Ar[r] * b[r]; });
+
+ return computed_ur;
+ }
+
+ NodeValuePerCell<const Rd>
+ _getFjr(const NodeValue<const Rd>& ur) const
+ {
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ const NodeValuePerCell<Rd> computed_Fjr(m_mesh.connectivity());
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ for (size_t r = 0; r < cell_nodes.size(); ++r) {
+ computed_Fjr(j, r) = m_Ajr(j, r) * (m_predicted_u[j] - ur[cell_nodes[r]]) + (m_predicted_p[j] * m_Djr(j, r));
+ }
+ });
- // std::cout << "Uk=" << Uk << '\n';
- // std::cout << "predicted_p" << m_predicted_p << '\n';
- // std::cout << m_predicted_u << '\n';
- // std::exit(0);
+ return computed_Fjr;
}
ImplicitAcousticSolver(const SolverType solver_type,
@@ -1647,6 +1851,11 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
return computed_Mj;
}();
+ m_u = u.cellValues();
+
+ m_predicted_u = m_u;
+ m_predicted_p = p.cellValues();
+
m_tau = [&]() {
CellValue<double> computed_tau(m_mesh.connectivity());
@@ -1673,7 +1882,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
if (m_is_implicit_cell[j]) {
- computed_g_1_exp_S_Cv_inv_g[j] = pow((gamma[j] - 1) * exp(entropy[j] / Cv[j]), m_inv_gamma[j]);
+ computed_g_1_exp_S_Cv_inv_g[j] = std::pow((gamma[j] - 1) * exp(entropy[j] / Cv[j]), m_inv_gamma[j]);
}
});
return computed_g_1_exp_S_Cv_inv_g;
@@ -1687,149 +1896,10 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
new_u = copy(u.cellValues());
new_E = copy(E.cellValues());
- _computeGradJU_AU(u, p, pi, dt);
-
- // u_j+1/2
- NodeValue<const Rd> ur = [&]() {
- const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
- const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
-
- NodeValue<Rd> b{m_mesh.connectivity()};
-
- parallel_for(
- m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
- const auto& node_to_cell = node_to_cell_matrix[r];
- const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
-
- Rd br = zero;
- for (size_t j = 0; j < node_to_cell.size(); ++j) {
- const CellId J = node_to_cell[j];
- const unsigned int R = node_local_number_in_its_cells[j];
- br += m_Ajr(J, R) * m_predicted_u[J] + m_predicted_p[J] * m_Djr(J, R);
- }
-
- b[r] = br;
- });
-
- for (const auto& boundary_condition : m_boundary_condition_list) {
- std::visit(
- [&](auto&& bc) {
- using T = std::decay_t<decltype(bc)>;
- if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
- if constexpr (Dimension == 1) {
- const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
-
- const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
- const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
-
- const auto& node_list = bc.faceList();
- const auto& value_list = bc.valueList();
- parallel_for(
- node_list.size(), PUGS_LAMBDA(size_t i_node) {
- const NodeId node_id = node_list[i_node];
- const auto& node_cell_list = node_to_cell_matrix[node_id];
- Assert(node_cell_list.size() == 1);
-
- CellId node_cell_id = node_cell_list[0];
- size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
-
- b[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
- });
- } else {
- const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
-
- const auto& face_to_cell_matrix = m_mesh.connectivity().faceToCellMatrix();
- const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
- const auto& face_local_numbers_in_their_cells = m_mesh.connectivity().faceLocalNumbersInTheirCells();
- const auto& face_cell_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
-
- const auto& face_list = bc.faceList();
- const auto& value_list = bc.valueList();
- for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
- const FaceId face_id = face_list[i_face];
- const auto& face_cell_list = face_to_cell_matrix[face_id];
- Assert(face_cell_list.size() == 1);
-
- CellId face_cell_id = face_cell_list[0];
- size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
-
- const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
-
- const auto& face_nodes = face_to_node_matrix[face_id];
-
- for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
- NodeId node_id = face_nodes[i_node];
- b[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
- }
- }
- }
- } else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
- // const auto cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
- auto is_boundary_node = m_mesh.connectivity().isBoundaryNode();
-
- const Rd& n = bc.outgoingNormal();
- const Rdxd I = identity;
- const Rdxd nxn = tensorProduct(n, n);
- const Rdxd P = I - nxn;
- const Array<const NodeId>& node_list = bc.nodeList();
- parallel_for(
- bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
- const NodeId r = node_list[r_number];
- b[r] = P * b[r];
- });
-
- } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- const auto& value_list = bc.valueList();
-
- parallel_for(
- node_list.size(), PUGS_LAMBDA(size_t i_node) {
- NodeId node_id = node_list[i_node];
- const auto& value = value_list[i_node];
- b[node_id] = value;
- });
- } else {
- throw UnexpectedError("boundary condition not handled 5");
- }
- },
- boundary_condition);
- }
-
- const NodeValue<Rd> computed_ur(m_mesh.connectivity());
- parallel_for(
- m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { computed_ur[r] = m_inv_Ar[r] * b[r]; });
-
- return computed_ur;
- }();
-
- // std::cout << "flux vitesse" << ur << '\n';
- // std::exit(0);
-
- // p_j+1/2
- const NodeValuePerCell<Rd>& Fjr = [&]() {
- const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
- const NodeValuePerCell<Rd> computed_Fjr(m_mesh.connectivity());
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
-
- for (size_t r = 0; r < cell_nodes.size(); ++r) {
- computed_Fjr(j, r) =
- m_Ajr(j, r) * (m_predicted_u[j] - ur[cell_nodes[r]]) + (m_predicted_p[j] * m_Djr(j, r));
- }
- });
-
- return computed_Fjr;
- }();
+ this->_computeGradJU_AU(u, p, pi, dt);
- // std::cout << "flux pression" << '\n';
- // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- // const size_t& nb_nodes = m_Ajr.numberOfSubValues(cell_id);
- // for (size_t r = 0; r < nb_nodes; ++r) {
- // std::cout << "Fjr(" << cell_id << "," << r << ")=" << Fjr(cell_id, r) << '\n';
- // }
- // }
+ NodeValue<const Rd> ur = this->_getUr();
+ NodeValuePerCell<const Rd> Fjr = this->_getFjr(ur);
// time n+1
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
--
GitLab