diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ce1a61e77e19542fadf80d6bc52d9618437af9b8..b950437fae7028efd644f555aca4e55ae1d16dae 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -621,6 +621,44 @@ SchemeModule::SchemeModule()
));
+ this
+ ->_addBuiltinFunction("implicit_eucclhyd_euler",
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& chi_explicit,
+ const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return ImplicitAcousticSolverHandler{ImplicitAcousticSolverHandler::SolverType::Eucclhyd,
+ rho,
+ c,
+ u,
+ p,
+ pi,
+ gamma,
+ Cv,
+ entropy,
+ bc_descriptor_list,
+ chi_explicit,
+ dt}
+ .solver()
+ .apply(dt, rho, u, E, c, p, pi, gamma, Cv, entropy);
+ }
+
+ ));
+
this->_addBuiltinFunction(
"implicit_eucclhyd_euler",
std::function(
@@ -655,6 +693,116 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("implicit_glace_euler",
+ 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 std::vector<std::shared_ptr<const IZoneDescriptor>>& explicit_zone_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return ImplicitAcousticSolverHandler{ImplicitAcousticSolverHandler::SolverType::
+ Glace2States,
+ rho,
+ c,
+ u,
+ p,
+ pi,
+ gamma,
+ Cv,
+ entropy,
+ bc_descriptor_list,
+ explicit_zone_list,
+ dt}
+ .solver()
+ .apply(dt, rho, u, E, c, p, pi, gamma, Cv, entropy);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("implicit_glace_euler",
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& chi_explicit,
+ const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return ImplicitAcousticSolverHandler{ImplicitAcousticSolverHandler::SolverType::
+ Glace2States,
+ rho,
+ c,
+ u,
+ p,
+ pi,
+ gamma,
+ Cv,
+ entropy,
+ bc_descriptor_list,
+ chi_explicit,
+ dt}
+ .solver()
+ .apply(dt, rho, u, E, c, p, pi, gamma, Cv, entropy);
+ }
+
+ ));
+
+ this->_addBuiltinFunction(
+ "implicit_glace_euler",
+ 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 IMesh>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> {
+ return ImplicitAcousticSolverHandler{ImplicitAcousticSolverHandler::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("acoustic_dt",
std::function(
diff --git a/src/scheme/ImplicitAcousticSolver.cpp b/src/scheme/ImplicitAcousticSolver.cpp
index d835bb7e101fc2236cd8dd255239afe85f354e3b..2b74689bf36b506d094ea05198707139a4007f3d 100644
--- a/src/scheme/ImplicitAcousticSolver.cpp
+++ b/src/scheme/ImplicitAcousticSolver.cpp
@@ -112,6 +112,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
NodeValuePerCell<const Rdxd> m_Ajr;
NodeValue<const Rdxd> m_inv_Ar;
+ NodeValuePerCell<const Rd> m_Djr;
+
CellValue<const Rd> m_predicted_u;
CellValue<const double> m_predicted_p;
@@ -152,15 +154,14 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
{
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
- const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
- const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+ const NodeValuePerCell<const Rd> Cjr_n = mesh_data.Cjr();
+ const NodeValuePerCell<const Rd> njr_n = mesh_data.njr();
NodeValuePerCell<Rdxd> Ajr{m_mesh.connectivity()};
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
switch (m_solver_type) {
case SolverType::Glace1State: {
- std::cout << "Glace 1 state" << '\n';
NodeValue<const double> rhocr = _getRhoCr(rho, c);
parallel_for(
@@ -168,25 +169,23 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const size_t& nb_nodes = Ajr.numberOfSubValues(j);
for (size_t r = 0; r < nb_nodes; ++r) {
const NodeId node_id = cell_to_node_matrix[j][r];
- Ajr(j, r) = tensorProduct(rhocr[node_id] * Cjr(j, r), njr(j, r));
+ Ajr(j, r) = tensorProduct(rhocr[node_id] * Cjr_n(j, r), njr_n(j, r));
}
});
break;
}
case SolverType::Glace2States: {
- std::cout << "Glace 2 states" << '\n';
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
const size_t& nb_nodes = Ajr.numberOfSubValues(j);
for (size_t r = 0; r < nb_nodes; ++r) {
- Ajr(j, r) = tensorProduct(rho[j] * c[j] * Cjr(j, r), njr(j, r));
+ Ajr(j, r) = tensorProduct(rho[j] * c[j] * Cjr_n(j, r), njr_n(j, r));
}
});
break;
}
case SolverType::Eucclhyd: {
- std::cout << "Eucclhyd" << '\n';
if constexpr (Dimension > 1) {
const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
@@ -224,7 +223,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
});
- // std::cout << "Nlr = " << Nlr << '\n';
break;
} else {
throw NotImplementedError("Eucclhyd switch is not implemented yet in 1d");
@@ -232,8 +230,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
- // std::cout << "Ajr = " << Ajr << '\n';
-
return Ajr;
}
@@ -246,17 +242,17 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
NodeValue<Rdxd> Ar{m_mesh.connectivity()};
parallel_for(
- m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
Rdxd sum = zero;
- 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);
+ const auto& node_to_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 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];
- sum += Ajr(J, R);
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const unsigned int i_node = node_local_number_in_its_cells[i_cell];
+ const CellId cell_id = node_to_cell[i_cell];
+ sum += Ajr(cell_id, i_node);
}
- Ar[r] = sum;
+ Ar[node_id] = sum;
});
NodeValue<int> boundary_node_s(m_mesh.connectivity());
@@ -306,6 +302,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
Ar[node_id] = identity;
});
+ } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
} else {
throw UnexpectedError("boundary condition not handled 1");
}
@@ -331,19 +328,19 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
int
mapP(CellId j) const
{
- // return (1 + Dimension) * m_implicit_cell_index[j];
- return m_implicit_cell_index[j];
+ return (1 + Dimension) * m_implicit_cell_index[j];
+ // return m_implicit_cell_index[j];
}
int
mapU(size_t k, CellId j) const
{
- // return (1 + Dimension) * m_implicit_cell_index[j] + k + 1;
- return m_number_of_implicit_cells + k + m_implicit_cell_index[j] * Dimension;
+ return (1 + Dimension) * m_implicit_cell_index[j] + k + 1;
+ // return m_number_of_implicit_cells + k + m_implicit_cell_index[j] * Dimension;
}
CRSMatrixDescriptor<double>
- _getA(NodeValuePerCell<const Rd> Cjr_half) const
+ _getA() const
{
count_getA++;
static bool is_get_A_started = false;
@@ -415,8 +412,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);
- const Rd coef =
- m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] * Cjr_half(id_cell_j, node_nb_in_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];
}
@@ -427,10 +423,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
- // std::cout << "sans cond limite=" << A << '\n';
-
- // NodeValue<Rdxd> inverse_ArxP{m_mesh.connectivity()};
-
NodeValue<int> boundary_node_s(m_mesh.connectivity());
boundary_node_s.fill(-1);
@@ -458,42 +450,40 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
// no changes for velocity conditions
- } // else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- // const auto& node_list = bc.faceList();
-
- // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
-
- // const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
-
- // const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
-
- // const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
-
- // for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- // for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- // const NodeId& node_id = node_list[i_node];
- // // if (m_is_implicit_node[node_id]) {
- // const auto& node_cell = node_to_cell_matrix[node_id];
- // 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 auto& node_local_number_in_its_cells =
- // node_local_numbers_in_their_cells.itemArray(node_id); const size_t node_nb_in_i =
- // node_local_number_in_its_cells[cell_i]; int index_p =
- // mapP(id_cell_i); int index_u = mapU(i_dimension, id_cell_i);
-
- // const Rd coef = m_Ajr(id_cell_i, node_local_number_in_its_cells[node_nb_in_i]) *
- // m_inv_Ar[node_id] *
- // Cjr(id_cell_i, node_local_number_in_its_cells[node_nb_in_i]);
- // A(index_p, index_u) += coef[i_dimension];
- // A(index_u, index_p) -= coef[i_dimension];
- // //}
- // }
- // }
- // }
- // }
- // }
- else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ // const auto& face_list = bc.faceList();
+
+ // const auto& face_local_numbers_in_their_cells = m_mesh.connectivity().faceLocalNumbersInTheirCells();
+
+ // const auto& face_to_cell_matrix = m_mesh.connectivity().faceToCellMatrix();
+
+ // for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ // for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ // const FaceId& face_id = face_list[i_face];
+ // // if (m_is_implicit_node[node_id]) {
+ // const auto& face_cell = face_to_cell_matrix[face_id];
+ // for (size_t cell_i = 0; cell_i < face_cell.size(); ++cell_i) {
+ // const CellId id_cell_i = face_cell[cell_i];
+ // if (m_is_implicit_cell[id_cell_i]) {
+ // // const auto& face_local_number_in_its_cells =
+ // // face_local_numbers_in_their_cells.itemArray(face_id); const size_t node_nb_in_i =
+ // // face_local_number_in_its_cells[cell_i];
+ // int index_p = mapP(id_cell_i);
+ // int index_u = mapU(i_dimension, id_cell_i);
+
+ // // const Rd coef = m_Ajr(id_cell_i, face_local_number_in_its_cells[node_nb_in_i]) *
+ // // m_inv_Ar[face_id] *
+ // // m_Djr(id_cell_i, face_local_number_in_its_cells[node_nb_in_i]);
+
+ // const Rd coef = zero;
+ // A(index_p, index_u) += coef[i_dimension];
+ // A(index_u, index_p) -= coef[i_dimension];
+ // //}
+ // }
+ // }
+ // }
+ // }
+ } else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
@@ -504,7 +494,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const Rdxd P = I - nxn;
NodeValue<Rdxd> inverse_ArxP{m_mesh.connectivity()};
- // std::cout << node_list << '\n';
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId& node_id = node_list[i_node];
if (boundary_node_s[node_id] != -1) {
@@ -526,7 +515,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const int j_index_p = mapP(id_cell_j);
const Rd coef =
- m_Ajr(id_cell_i, node_nb_in_i) * inverse_ArxP[node_id] * Cjr_half(id_cell_j, node_nb_in_j);
+ m_Ajr(id_cell_i, node_nb_in_i) * inverse_ArxP[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];
}
@@ -544,8 +533,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
get_A_t.pause();
- // std::cout << "avec cond limite" << A << '\n';
-
return A;
}
@@ -600,12 +587,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
computed_tau_iter[j] = m_g_1_exp_S_Cv_inv_g[j] * pow(-Uk[k] + pi[j], -m_inv_gamma[j]);
}
});
- // for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
- // // std::cout << "exp=" << m_g_1_exp_S_Cv_inv_g[j] << '\n';
- // // std::cout << "pi=" << pi[j] << '\n';
- // std::cout << "tau_iter=" << computed_tau_iter[j] << '\n';
- // std::cout << "Un=" << Un << '\n';
- // }
return computed_tau_iter;
}();
@@ -615,16 +596,13 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
computed_gradJ = zero;
const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
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();
auto is_boundary_node = m_mesh.connectivity().isBoundaryNode();
for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
if (not is_boundary_node[node_id]) {
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);
- Rd sum_Cjrpj = zero;
+ Rd sum_Djrpj = zero;
Rd sum_Ajruj = zero;
Rd Uk_i = zero;
Rd Uk_j = zero;
@@ -633,7 +611,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_cell[id_cell_i]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
const size_t i_index_p = mapP(id_cell_i);
- sum_Cjrpj += Uk[i_index_p] * Cjr(id_cell_i, node_nb_in_i);
+ sum_Djrpj += Uk[i_index_p] * m_Djr(id_cell_i, node_nb_in_i);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const size_t i_index_u = mapU(i_dimension, id_cell_i);
@@ -649,7 +627,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const size_t node_nb_in_j = node_local_number_in_its_cells[cell_j];
const size_t j_index_p = mapP(id_cell_j);
- computed_gradJ[j_index_p] += dot(Cjr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id] * sum_Cjrpj);
+ computed_gradJ[j_index_p] += dot(m_Djr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id] * sum_Djrpj);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const size_t j_index_u = mapU(i_dimension, id_cell_j);
@@ -682,12 +660,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_node[node_id]) {
if (not is_boundary_node[node_id]) {
const auto& node_cells = node_to_cell_matrix[node_id];
- Rd sum_Cjrpj_exp = zero;
+ Rd sum_Djrpj_exp = zero;
Rd sum_Ajruj_exp = zero;
for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
if (not m_is_implicit_cell[node_cells[i_cell]]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[i_cell];
- sum_Cjrpj_exp += p[node_cells[i_cell]] * Cjr(node_cells[i_cell], node_nb_in_i);
+ sum_Djrpj_exp += p[node_cells[i_cell]] * m_Djr(node_cells[i_cell], node_nb_in_i);
sum_Ajruj_exp += m_Ajr(node_cells[i_cell], node_nb_in_i) * u[node_cells[i_cell]];
}
}
@@ -696,13 +674,13 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_cell[id_cell_j]) {
const size_t node_nb_in_j = node_local_number_in_its_cells[j_cell];
const size_t j_index_p = mapP(id_cell_j);
- computed_gradJ[j_index_p] -= dot(Cjr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id] * sum_Cjrpj_exp) +
- dot(Cjr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id] * sum_Ajruj_exp);
+ computed_gradJ[j_index_p] -= dot(m_Djr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id] * sum_Djrpj_exp) +
+ dot(m_Djr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id] * sum_Ajruj_exp);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const size_t j_index_u = mapU(i_dimension, id_cell_j);
computed_gradJ[j_index_u] -=
- (m_Ajr(id_cell_j, node_nb_in_j) * m_inv_Ar[node_id] * sum_Cjrpj_exp +
+ (m_Ajr(id_cell_j, node_nb_in_j) * m_inv_Ar[node_id] * sum_Djrpj_exp +
m_Ajr(id_cell_j, node_nb_in_j) * m_inv_Ar[node_id] * sum_Ajruj_exp)[i_dimension];
}
}
@@ -710,9 +688,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
}
- // std::cout << "gradJ sans cond limit=" << computed_gradJ << '\n';
- // NodeValue<int> boundary_node_sym(m_mesh.connectivity());
- // boundary_node_sym.fill(-1);
NodeValue<int> boundary_node_s(m_mesh.connectivity());
boundary_node_s.fill(-1);
@@ -745,10 +720,9 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// const auto& value_list = bc.valueList();
// MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
- // const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
// const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
// const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
- // Rd sum_Cjrpext = zero;
+ // Rd sum_Djrpext = zero;
// for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
// const NodeId& node_id = node_list[i_node];
@@ -761,7 +735,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// node_local_numbers_in_their_cells.itemArray(node_id); const size_t node_nb_in_i =
// node_local_number_in_its_cells[cell_i];
- // sum_Cjrpext += value_list[i_node] * Cjr(id_cell_i, node_nb_in_i);
+ // sum_Djrpext += value_list[i_node] * m_Djr(id_cell_i, node_nb_in_i);
// }
// }
// for (size_t cell_j = 0; cell_j < node_cell.size(); ++cell_j) {
@@ -770,11 +744,11 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// const auto& node_local_number_in_its_cells =
// node_local_numbers_in_their_cells.itemArray(node_id); const size_t node_nb_in_j =
// node_local_number_in_its_cells[cell_j]; int index_p =
- // mapP(id_cell_j); computed_gradJ[index_p] += dot(Cjr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id]
- // * sum_Cjrpext); for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ // mapP(id_cell_j); computed_gradJ[index_p] += dot(m_Djr(id_cell_j, node_nb_in_j), m_inv_Ar[node_id]
+ // * sum_Djrpext); for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
// int index_u = mapU(i_dimension, id_cell_j);
// computed_gradJ[index_u] +=
- // (m_Ajr(id_cell_j, node_nb_in_j) * m_inv_Ar[node_id] * sum_Cjrpext)[i_dimension];
+ // (m_Ajr(id_cell_j, node_nb_in_j) * m_inv_Ar[node_id] * sum_Djrpext)[i_dimension];
// }
// }
// }
@@ -782,8 +756,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// }
if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
- const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
const Rd& n = bc.outgoingNormal();
@@ -800,9 +772,9 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const auto& node_cell = node_to_cell_matrix[node_id];
inverse_ArxP[node_id] = m_inv_Ar[node_id] * P;
- Rd sum_Cjrpj = zero;
+ Rd sum_Djrpj = zero;
Rd sum_Ajruj = zero;
- Rd sum_Cjrpj_exp = zero;
+ Rd sum_Djrpj_exp = zero;
Rd sum_Ajruj_exp = zero;
Rd Uk_i = zero;
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
@@ -812,7 +784,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_is_implicit_cell[id_cell_i]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
const size_t i_index_p = mapP(id_cell_i);
- sum_Cjrpj += Uk[i_index_p] * Cjr(id_cell_i, node_nb_in_i);
+ sum_Djrpj += Uk[i_index_p] * m_Djr(id_cell_i, node_nb_in_i);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const size_t i_index_u = mapU(i_dimension, id_cell_i);
@@ -825,7 +797,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
for (size_t i_cell = 0; i_cell < node_cell.size(); ++i_cell) {
if (not m_is_implicit_cell[node_cell[i_cell]]) {
const size_t node_nb_in_i = node_local_number_in_its_cells[i_cell];
- sum_Cjrpj_exp += p[node_cell[i_cell]] * Cjr(node_cell[i_cell], node_nb_in_i);
+ sum_Djrpj_exp += p[node_cell[i_cell]] * m_Djr(node_cell[i_cell], node_nb_in_i);
sum_Ajruj_exp += m_Ajr(node_cell[i_cell], node_nb_in_i) * u[node_cell[i_cell]];
}
}
@@ -837,9 +809,9 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
int index_p = mapP(id_cell_j);
computed_gradJ[index_p] +=
- dot(Cjr(id_cell_j, node_nb_in_j), inverse_ArxP[node_id] * sum_Cjrpj) -
- dot(Cjr(id_cell_j, node_nb_in_j), inverse_ArxP[node_id] * sum_Cjrpj_exp) -
- dot(Cjr(id_cell_j, node_nb_in_j), inverse_ArxP[node_id] * sum_Ajruj_exp);
+ dot(m_Djr(id_cell_j, node_nb_in_j), inverse_ArxP[node_id] * sum_Djrpj) -
+ dot(m_Djr(id_cell_j, node_nb_in_j), inverse_ArxP[node_id] * sum_Djrpj_exp) -
+ dot(m_Djr(id_cell_j, node_nb_in_j), inverse_ArxP[node_id] * sum_Ajruj_exp);
;
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
int j_index_u = mapU(i_dimension, id_cell_j);
@@ -851,7 +823,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
computed_gradJ[index_u] +=
(-m_Ajr(id_cell_j, node_nb_in_j) * inverse_ArxP[node_id] * sum_Ajruj +
m_Ajr(id_cell_j, node_nb_in_j) * Uk_j)[i_dimension] -
- (m_Ajr(id_cell_j, node_nb_in_j) * inverse_ArxP[node_id] * sum_Cjrpj_exp +
+ (m_Ajr(id_cell_j, node_nb_in_j) * inverse_ArxP[node_id] * sum_Djrpj_exp +
m_Ajr(id_cell_j, node_nb_in_j) * inverse_ArxP[node_id] * sum_Ajruj_exp)[i_dimension];
}
}
@@ -862,8 +834,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const auto& node_list = bc.nodeList();
const auto& value_list = bc.valueList();
- MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
- const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
@@ -881,9 +851,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
int i_index_u = mapU(i_dimension, id_cell_i);
Uk_i[i_dimension] = Uk[i_index_u];
- // std::cout << "Uk_i[" << i_dimension << "]=" << Uk_i[i_dimension] << '\n';
- // std::cout << "Uk=" << Uk << '\n';
- // std::cout << "Uk[" << i_index_u << "]=" << Uk[i_index_u] << '\n';
}
}
}
@@ -894,32 +861,23 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
int i_index_p = mapP(id_cell_i);
- computed_gradJ[i_index_p] += -dot(Cjr(id_cell_i, node_nb_in_i), value_list[i_node]);
- // std::cout << "gradJp[" << i_index_p << "]=" << computed_gradJ[i_index_p] << '\n';
- // std::cout << "Cjr(" << id_cell_i << "," << node_nb_in_i << ")=" << Cjr(id_cell_i, node_nb_in_i)
- // << '\n'; std::cout << "value[" << i_node << "]=" << value_list[i_node] << '\n';
+ computed_gradJ[i_index_p] += -dot(m_Djr(id_cell_i, node_nb_in_i), value_list[i_node]);
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
int i_index_u = mapU(i_dimension, id_cell_i);
- // Uk_i[i_dimension] = Uk[i_index_u];
- // std::cout << "Uk_i" << Uk_i << '\n';
computed_gradJ[i_index_u] +=
(m_Ajr(id_cell_i, node_nb_in_i) * (Uk_i - value_list[i_node]))[i_dimension];
- // std::cout << "Ajr(" << id_cell_i << "," << node_nb_in_i << ")=" << m_Ajr(id_cell_i,
- // node_nb_in_i) << '\n'; std::cout << "uk-value[" << i_node << "]=" << Uk_i - value_list[i_node]
- // << '\n'; std::cout << "gradJu[" << i_index_u << "]=" << computed_gradJ[i_index_u] << '\n';
}
}
}
}
}
+ } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
} else {
throw UnexpectedError("boundary condition not handled 3");
}
},
boundary_condition);
}
- // std::cout << "gradJ=" << computed_gradJ << '\n';
- // std::exit(0);
return computed_gradJ;
}
@@ -934,13 +892,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
{
Vector<double> gradJ = this->_getGradJ(Un, Uk, p, pi, u, dt);
- // std::cout << "gradJ" << gradJ << '\n';
- // CRSMatrix A_crs{A.getCRSMatrix()};
- // std::exit(0);
- // Vector<double> AU = A_crs * Uk;
Vector<double> AU = A * Uk;
- // std::cout << "AU" << AU << '\n';
Vector<double> F = gradJ - AU;
return F;
@@ -990,7 +943,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
nb_neighboors.clear();
}
}
- // std::cout << non_zeros << '\n';
CRSMatrixDescriptor Hess_J{(Dimension + 1) * m_number_of_implicit_cells,
(Dimension + 1) * m_number_of_implicit_cells, non_zeros};
@@ -1028,58 +980,58 @@ 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 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_p = mapP(id_cell_i);
+ for (size_t j_cell = 0; j_cell < node_cell.size(); ++j_cell) {
+ const CellId cell_id_j = node_cell[j_cell];
+ if (m_is_implicit_cell[cell_id_j]) {
+ 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(id_cell_i, node_nb_in_i);
+ const auto invArCir = inv_Ar * Cjr(cell_id_j, node_nb_in_j_cell);
- 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);
- Hess_J(i_index_p, j_index_p) += dot(invArCir, Cjr(id_cell_j, node_nb_in_j));
+ 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));
}
}
}
}
- 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 j_cell = 0; j_cell < node_cell.size(); ++j_cell) {
+ const CellId cell_id_j = node_cell[j_cell];
+ if (m_is_implicit_cell[cell_id_j]) {
+ const size_t node_nb_in_j_cell = node_local_number_in_its_cells[j_cell];
- const auto Ajr = m_Ajr(id_cell_i, node_nb_in_i);
+ const auto Ajr = m_Ajr(cell_id_j, node_nb_in_j_cell);
- for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- const int i_index_u = mapU(i_dimension, id_cell_i);
- Hess_J(i_index_u, i_index_u) += Ajr(i_dimension, i_dimension);
+ for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
+ const int j_index_u = mapU(j_dimension, cell_id_j);
+ Hess_J(j_index_u, j_index_u) += Ajr(j_dimension, j_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 j_cell = 0; j_cell < node_cell.size(); ++j_cell) {
+ const CellId cell_id_j = node_cell[j_cell];
+ if (m_is_implicit_cell[cell_id_j]) {
+ const size_t node_nb_in_j_cell = node_local_number_in_its_cells[j_cell];
- const auto Air_invAr = m_Ajr(id_cell_i, node_nb_in_i) * inv_Ar;
+ const auto Ajr_invAr = m_Ajr(cell_id_j, node_nb_in_j_cell) * inv_Ar;
- 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];
+ 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_i_cell = node_local_number_in_its_cells[i_cell];
- const auto Air_invAr_Ajr = Air_invAr * m_Ajr(id_cell_j, node_nb_in_j);
+ const auto Ajr_invAr_Air = Ajr_invAr * m_Ajr(cell_id_i, node_nb_in_i_cell);
- for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- const int i_index_u = mapU(i_dimension, id_cell_i);
+ for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
+ const int index_u_j = mapU(j_dimension, cell_id_j);
for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
- const int j_index_u = mapU(j_dimension, id_cell_j);
- Hess_J(i_index_u, j_index_u) -= Air_invAr_Ajr(i_dimension, j_dimension);
+ const int index_u_i = mapU(j_dimension, cell_id_i);
+ Hess_J(index_u_j, index_u_i) -= Ajr_invAr_Air(j_dimension, j_dimension);
}
}
}
@@ -1224,70 +1176,68 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
}
- }
- // else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- // const auto& node_list = bc.faceList();
-
- // const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
- // const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
-
- // for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- // const NodeId& node_id = node_list[i_node];
- // // if (m_is_implicit_node[node_id]) {
- // const auto& node_cell = node_to_cell_matrix[node_id];
-
- // 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 auto& node_local_number_in_its_cells =
- // node_local_numbers_in_their_cells.itemArray(node_id);
-
- // const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
- // const int i_index_p = mapP(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);
-
- // 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));
- // }
- // }
- // }
- // }
- // 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 auto& node_local_number_in_its_cells =
- // node_local_numbers_in_their_cells.itemArray(node_id);
-
- // const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
- // for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- // const int i_index_u = mapU(i_dimension, id_cell_i);
-
- // Hess_J(i_index_u, i_index_u) += m_Ajr(id_cell_i, node_nb_in_i)(i_dimension, i_dimension);
-
- // 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];
- // for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
- // const int j_index_u = mapU(j_dimension, id_cell_j);
-
- // Hess_J(i_index_u, j_index_u) += -(m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] *
- // m_Ajr(id_cell_j, node_nb_in_j))(i_dimension,
- // j_dimension);
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- // }
- else {
+ } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ // const auto& node_list = bc.faceList();
+
+ // const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
+ // const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+ // for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ // const NodeId& node_id = node_list[i_node];
+ // // if (m_is_implicit_node[node_id]) {
+ // const auto& node_cell = node_to_cell_matrix[node_id];
+
+ // 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 auto& node_local_number_in_its_cells =
+ // node_local_numbers_in_their_cells.itemArray(node_id);
+
+ // const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
+ // const int i_index_p = mapP(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);
+
+ // 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));
+ // }
+ // }
+ // }
+ // }
+ // 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 auto& node_local_number_in_its_cells =
+ // node_local_numbers_in_their_cells.itemArray(node_id);
+
+ // const size_t node_nb_in_i = node_local_number_in_its_cells[cell_i];
+ // for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ // const int i_index_u = mapU(i_dimension, id_cell_i);
+
+ // Hess_J(i_index_u, i_index_u) += m_Ajr(id_cell_i, node_nb_in_i)(i_dimension, i_dimension);
+
+ // 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];
+ // for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
+ // const int j_index_u = mapU(j_dimension, id_cell_j);
+
+ // Hess_J(i_index_u, j_index_u) += -(m_Ajr(id_cell_i, node_nb_in_i) * m_inv_Ar[node_id] *
+ // m_Ajr(id_cell_j, node_nb_in_j))(i_dimension,
+ // j_dimension);
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ } else {
throw UnexpectedError("boundary condition not handled 4");
}
},
@@ -1312,10 +1262,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CRSMatrixDescriptor<double> Hess_J = this->_getHessianJ(dt, pi, Uk);
- // std::cout << "hessJ" << Hess_J_crs << '\n';
- // std::exit(0);
- // CRSMatrix A_crs{A.getCRSMatrix()};
-
CRSMatrix Hess_J_crs = Hess_J.getCRSMatrix();
CRSMatrix gradient_f = Hess_J_crs - A;
@@ -1406,8 +1352,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
_computeGradJU_AU(const DiscreteVectorFunction& u,
const DiscreteScalarFunction& p,
const DiscreteScalarFunction& pi,
- const double& dt,
- const NodeValuePerCell<const Rd>& Cjr_half)
+ const double& dt)
{
// const auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
const auto cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -1429,12 +1374,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
return implicit_cell_index;
}();
- // std::cout << "cellules implicites=true" << '\n';
- // std::cout << std::boolalpha << m_is_implicit_cell << '\n';
- // std::cout << "renumerotation cellules implicites" << '\n';
- // std::cout << m_implicit_cell_index << '\n';
- // std::cout << "nb implicit cells" << m_number_of_implicit_cells << '\n';
-
m_is_implicit_node = [&] {
NodeValue<bool> is_implicit_node(m_mesh.connectivity());
is_implicit_node.fill(false);
@@ -1468,25 +1407,13 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
return implicit_node_index;
}();
- // std::cout << "noeuds implicites=true" << '\n';
- // std::cout << std::boolalpha << m_is_implicit_node << '\n';
- // std::cout << "renumerotation noeuds implicites" << '\n';
- // std::cout << m_implicit_node_index << '\n';
-
- // std::cout << "inv_Ar" << '\n';
- // std::cout << m_inv_Ar << '\n';
-
- // const CRSMatrixDescriptor<double>& A = this->_getA(Cjr_half);
+ std::cout << "building A: ";
+ const CRSMatrix A = this->_getA().getCRSMatrix();
+ std::cout << "done\n";
- const CRSMatrix A = this->_getA(Cjr_half).getCRSMatrix();
-
- // std::cout << "A" << A << '\n';
- // std::exit(0);
Vector<double> Un = this->_getU(p, u);
Vector<double> Uk = copy(Un);
- // std::cout << "Un:" << Un << '\n';
-
int nb_iter = 0;
double norm_inf_sol;
@@ -1500,25 +1427,44 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
double norm_inf_Un = max(abs_Un);
if (m_number_of_implicit_cells > 0) {
do {
- // std::cout << "iteration=" << nb_iter << '\n';
nb_iter++;
- std::cout << "nb_iter = " << nb_iter << '\n';
Vector<double> f = this->_getF(A, Un, Uk, p, pi, u, dt);
- // std::exit(0);
+
+ std::cout << "building gradf: ";
CRSMatrix<double> gradient_f = this->_getGradientF(A, Uk, pi, dt);
+ std::cout << " done\n";
+
+ // std::cout << "grad f = " << gradient_f << '\n';
+ // std::cout << "f=" << f << '\n';
+ std::cout << " -> f: [" << min(f) << ", " << max(f) << "]\n";
+ // std::cout << "grad f=" << gradient_f << '\n';
+
+ auto l2Norm = [](const Vector<double>& x) {
+ double sum2 = 0;
+ for (size_t i = 0; i < x.size(); ++i) {
+ sum2 += x[i] * x[i];
+ }
+ return std::sqrt(sum2);
+ };
- Vector<double> sol = copy(Uk);
+ Vector<double> sol{Uk.size()};
+ sol.fill(0);
static bool solver_is_init = false;
if (not solver_is_init) {
solver_t.stop();
solver_is_init = true;
}
+
solver_t.start();
LinearSolver solver;
solver.solveLocalSystem(gradient_f, sol, f);
solver_t.pause();
+ std::cout << " -> sol: [" << min(sol) << ", " << max(sol) << "]\n";
+ std::cout << "l2Norm(f) = " << l2Norm(f) << '\n';
+ std::cout << "l2Norm(sol) = " << l2Norm(sol) << '\n';
+ 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;
@@ -1533,28 +1479,32 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
Uk = U_next;
- std::cout << "ratio Newton = " << norm_inf_sol / norm_inf_Un << "\n";
+ std::cout << "nb_iter= " << nb_iter << " | ratio Newton = " << norm_inf_sol / norm_inf_Un << "\n";
// when it is a hard case we relax newton
- // bool neg_pressure = false;
- // for (CellId j = 0; j < m_number_of_implicit_cells; ++j) {
- // size_t k = mapP(j);
- // if (-U_next[k] + pi[j] < 0) {
- // U_next[k] = Uk[k];
- // // std::cout << "p est negatif pour la maille " << j << " -Uk[k] + pi[j]=" << -U_next[k] + pi[j] <<
- // '\n'; neg_pressure = true;
- // }
- // }
-
+ 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) {
+ ++neg_pressure_count;
+ min_pressure = std::min(min_pressure, -U_next[k] + pi[j]);
+ U_next[k] = Uk[k] + sol[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';
+ }
// Uk -= sol;
// if (neg_pressure) {
// std::cout << "HAS NEG PRESSURE, RELAXING\n";
- // Uk -= 0.01 * sol;
- // // } else {
- // // Uk = U_next;
+ // Uk -= 0.01 * sol;
+ // } else {
+ // Uk -= sol;
// }
- } while ((norm_inf_sol > 1e-6 * norm_inf_Un) and (nb_iter < 200));
+ } while ((norm_inf_sol > 1e-15 * norm_inf_Un) and (nb_iter < 200));
}
m_predicted_u = [&] {
CellValue<Rd> predicted_u = copy(u.cellValues());
@@ -1629,6 +1579,34 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}();
}
+ ImplicitAcousticSolver(const SolverType solver_type,
+ const std::shared_ptr<const MeshType>& p_mesh,
+ const DiscreteScalarFunction&,
+ const DiscreteScalarFunction&,
+ const DiscreteVectorFunction&,
+ const DiscreteScalarFunction&,
+ const DiscreteScalarFunction&,
+ const DiscreteScalarFunction&,
+ const DiscreteScalarFunction&,
+ const DiscreteScalarFunction&,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const DiscreteScalarFunction& chi_explicit,
+ const double&)
+ : m_solver_type{solver_type},
+ m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)},
+ m_mesh{*p_mesh}
+ {
+ m_is_implicit_cell = [&] {
+ CellValue<bool> is_implicit_cell(m_mesh.connectivity());
+
+ parallel_for(
+ m_mesh.numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { is_implicit_cell[cell_id] = (chi_explicit[cell_id] == 0); });
+
+ return is_implicit_cell;
+ }();
+ }
+
public:
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -1653,7 +1631,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
const NodeValuePerCell<const Rd> Cjr_n = mesh_data.Cjr();
- NodeValuePerCell<Rd> Cjr_half = copy(Cjr_n);
+ m_Djr = copy(Cjr_n);
CellValue<double> new_rho = copy(rho.cellValues());
CellValue<Rd> new_u = copy(u.cellValues());
@@ -1679,18 +1657,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
m_Ajr = this->_computeAjr(rho, c);
- // 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 << "Ajr(" << cell_id << "," << r << ")=" << m_Ajr(cell_id, r) << '\n';
- // }
- // }
-
NodeValue<const Rdxd> Ar = this->_computeAr(m_Ajr);
- // std::cout << "Ar=" << '\n';
- // std::cout << Ar << '\n';
-
m_inv_Ar = this->_computeInvAr(Ar);
m_inv_gamma = [&]() {
@@ -1719,7 +1687,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
new_u = copy(u.cellValues());
new_E = copy(E.cellValues());
- _computeGradJU_AU(u, p, pi, dt, Cjr_half);
+ _computeGradJU_AU(u, p, pi, dt);
// u_j+1/2
NodeValue<const Rd> ur = [&]() {
@@ -1737,7 +1705,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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] * Cjr_n(J, R);
+ br += m_Ajr(J, R) * m_predicted_u[J] + m_predicted_p[J] * m_Djr(J, R);
}
b[r] = br;
@@ -1848,7 +1816,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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] * Cjr_n(j, r));
+ m_Ajr(j, r) * (m_predicted_u[j] - ur[cell_nodes[r]]) + (m_predicted_p[j] * m_Djr(j, r));
}
});
@@ -1889,24 +1857,25 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
new_u[j] -= dt_over_Mj * momentum_fluxes;
new_E[j] -= dt_over_Mj * energy_fluxes;
});
- // std::cout << "newton converge, mise a jour" << '\n';
- // std::cout << "rho=" << new_rho << '\n';
- // std::cout << MeshDataManager::instance().getMeshData(*mesh).xj() << '\n';
- // std::cout << nouveau_Vj << '\n';
-
- CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
- // std::cout << "volume mis a jour" << '\n';
- // update Cjr_half
+ // update Djr
const NodeValuePerCell<const Rd> new_Cjr = MeshDataManager::instance().getMeshData(*new_mesh).Cjr();
+ CellValue<double> new_Vj{m_mesh.connectivity()};
parallel_for(
- mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+ mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ new_Vj[j] = 0;
+ auto cell_nodes = cell_to_node_matrix[j];
- for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
- // std::cout << "rho(" << j << ")=" << rho[j] << '\n';
- // std::cout << "new_rho(" << j << ")=" << new_rho[j] << '\n';
- }
+ 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]; });
CellValue<double> new_tau(m_mesh.connectivity());
@@ -1915,7 +1884,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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(Cjr_half(j, R), ur[r]);
+ new_tau_fluxes += dot(m_Djr(j, R), ur[r]);
}
new_tau[j] = m_tau[j] + (dt / m_Mj[j]) * new_tau_fluxes;
// std::cout << "new_tau(" << j << ")=" << new_tau[j] << '\n';
@@ -1928,33 +1897,131 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// 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]);
- max_tau_error = std::max(max_tau_error, std::abs(new_tau[j] - 1. / new_rho[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 << "max_tau_error=" << max_tau_error << '\n';
- for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
- 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];
- Cjr_half(j, R) = 0.5 * (Cjr_n(j, R) + new_Cjr(j, R));
- // Cjr_half(j, R) = Cjr_n(j, R);
- // std::cout << "Cjr_n(" << j << "," << R << ")=" << Cjr_n(j, R) << '\n';
- // std::cout << "Cjr_new(" << j << "," << R << ")=" << new_Cjr(j, R) << '\n';
- // std::cout << "Cjr_half(" << j << "," << R << ")=" << Cjr_half(j, R) << '\n';
- }
+ 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);
+ }
+ }
+ });
+
+ 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;
+
+ 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) {
+ const Rd two_dxs = 2 * dxr_n[s];
+ Nr +=
+ crossProduct(2 * ((dxr_np1[r] - 2 * dxr_np1[s]) - (dxr_n[r] - 2 * dxr_n[s])), xr_n[face_nodes[s]]);
+ Nr -= crossProduct((1. / 6) * ((2 * dxr_n[r] - dxr_n[s]) - (2 * dxr_n[r] - 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_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ 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);
+
+ 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);
+
+ 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;
+ }
+ }
+ 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);
+ }
+ }
+ }
+ });
+
+ 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);
+ }
+ }
+ });
+
+ m_Djr = Djr;
+ }
+
+ if constexpr (Dimension > 1) {
+ std::cout << "number_iter_Djr=" << number_iter << " max_tau_error=" << max_tau_error << '\n';
}
- std::cout << "number_iter_Cjr=" << number_iter << '\n';
// std::cout << "new rho=" << new_rho << '\n';
- } while ((max_tau_error > 1e-13) and (number_iter < 2));
+ } 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 << ")"
- << " getF=" << getF_t << "(" << count_getF << ")"
- << " getGradF=" << getGradF_t << "(" << count_getGradF << ")"
- << " HJ_A=" << HJ_A_t << '\n';
- std::cout << "solver=" << solver_t << " total= " << implicit_t << '\n';
+ // std::cout << "getA=" << get_A_t << "(" << count_getA << ")"
+ // << " getF=" << getF_t << "(" << count_getF << ")"
+ // << " getGradF=" << getGradF_t << "(" << count_getGradF << ")"
+ // << " HJ_A=" << HJ_A_t << '\n';
+ // std::cout << "solver=" << solver_t << " total= " << implicit_t << '\n';
return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction{new_mesh, new_rho}),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction{new_mesh, new_u}),
@@ -2020,6 +2087,34 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
dt)
{}
+ ImplicitAcousticSolver(const SolverType solver_type,
+ const std::shared_ptr<const IMesh>& mesh,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& chi_explicit,
+ const double& dt)
+ : ImplicitAcousticSolver(solver_type,
+ std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
+ rho->get<DiscreteScalarFunction>(),
+ c->get<DiscreteScalarFunction>(),
+ u->get<DiscreteFunctionP0<Dimension, const Rd>>(),
+ p->get<DiscreteScalarFunction>(),
+ pi->get<DiscreteScalarFunction>(),
+ gamma->get<DiscreteScalarFunction>(),
+ Cv->get<DiscreteScalarFunction>(),
+ entropy->get<DiscreteScalarFunction>(),
+ bc_descriptor_list,
+ chi_explicit->get<DiscreteScalarFunction>(),
+ dt)
+ {}
+
ImplicitAcousticSolver() = default;
ImplicitAcousticSolver(ImplicitAcousticSolver&&) = default;
~ImplicitAcousticSolver() = default;
@@ -2190,3 +2285,51 @@ ImplicitAcousticSolverHandler::ImplicitAcousticSolverHandler(
}
}
}
+
+ImplicitAcousticSolverHandler::ImplicitAcousticSolverHandler(
+ const SolverType solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& chi_explicit,
+ const double& dt)
+{
+ std::shared_ptr i_mesh = getCommonMesh({rho, c, u, p, pi, gamma, Cv, chi_explicit});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho, c, u, p}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ switch (i_mesh->dimension()) {
+ case 1: {
+ m_implicit_acoustic_solver =
+ std::make_unique<ImplicitAcousticSolver<1>>(solver_type, i_mesh, rho, c, u, p, pi, gamma, Cv, entropy,
+ bc_descriptor_list, chi_explicit, dt);
+ break;
+ }
+ case 2: {
+ m_implicit_acoustic_solver =
+ std::make_unique<ImplicitAcousticSolver<2>>(solver_type, i_mesh, rho, c, u, p, pi, gamma, Cv, entropy,
+ bc_descriptor_list, chi_explicit, dt);
+ break;
+ }
+ case 3: {
+ m_implicit_acoustic_solver =
+ std::make_unique<ImplicitAcousticSolver<3>>(solver_type, i_mesh, rho, c, u, p, pi, gamma, Cv, entropy,
+ bc_descriptor_list, chi_explicit, dt);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ }
+}
diff --git a/src/scheme/ImplicitAcousticSolver.hpp b/src/scheme/ImplicitAcousticSolver.hpp
index 2ad0e9f3a66f5f172ff78f808356ea8da69b5f18..c3c7eb4feccf9cedc1ede2fc239ae1ea6805dfdf 100644
--- a/src/scheme/ImplicitAcousticSolver.hpp
+++ b/src/scheme/ImplicitAcousticSolver.hpp
@@ -90,6 +90,20 @@ class ImplicitAcousticSolverHandler
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const std::vector<std::shared_ptr<const IZoneDescriptor>>& explicit_zone_list,
const double& dt);
+
+ ImplicitAcousticSolverHandler(
+ const SolverType solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& chi_explicit,
+ const double& dt);
};
#endif // IMPLICIT_ACOUSTIC_SOLVER_HPP