diff --git a/src/scheme/ImplicitAcousticSolver.cpp b/src/scheme/ImplicitAcousticSolver.cpp
index 7cfed87083e24a82ecdac253669fd58fa1bf382e..411c4112c0f1af52061ffb791094732cdd03f538 100644
--- a/src/scheme/ImplicitAcousticSolver.cpp
+++ b/src/scheme/ImplicitAcousticSolver.cpp
@@ -5,7 +5,9 @@
#include <mesh/MeshCellZone.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshLineNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
@@ -13,20 +15,28 @@
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <output/NamedItemValueVariant.hpp>
+#include <output/VTKWriter.hpp>
#include <utils/Timer.hpp>
#include <variant>
#include <vector>
+// VTKWriter vtk_writer("imp/debug", 0);
+
+int count_newton = 0;
+int count_Djr = 0;
+int count_timestep = 0;
+
Timer solver_t;
Timer get_A_t;
Timer getF_t;
Timer getGradF_t;
Timer HJ_A_t;
-static int count_getA = 0;
-static int count_getF = 0;
-static int count_getGradF = 0;
+int count_getA = 0;
+int count_getGradF = 0;
template <size_t Dimension>
double
@@ -89,12 +99,13 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;
+ class AxisBoundaryCondition;
class PressureBoundaryCondition;
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
- using BoundaryCondition =
- std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+ using BoundaryCondition = std::
+ variant<AxisBoundaryCondition, PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -183,7 +194,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
Ajr(j, r) = tensorProduct(rho[j] * c[j] * Cjr_n(j, r), njr_n(j, r));
}
});
-
break;
}
case SolverType::Eucclhyd: {
@@ -256,24 +266,61 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
Ar[node_id] = sum;
});
- NodeValue<int> boundary_node_s(m_mesh.connectivity());
- boundary_node_s.fill(-1);
+ NodeValue<bool> has_boundary_condition{m_mesh.connectivity()};
+ has_boundary_condition.fill(false);
+ // velocity bc
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<VelocityBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId& node_id = node_list[i_node];
- boundary_node_s[node_id] = 1;
- }
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ if (not has_boundary_condition[node_id]) {
+ Ar[node_id] = identity;
+
+ has_boundary_condition[node_id] = true;
+ }
+ });
}
},
boundary_condition);
}
+ if constexpr (Dimension > 1) {
+ // axis bc
+ 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<AxisBoundaryCondition, T>) {
+ const Rd& t = bc.direction();
+
+ const Rdxd I = identity;
+ const Rdxd txt = tensorProduct(t, t);
+
+ const Array<const NodeId>& node_list = bc.nodeList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId& node_id = node_list[i_node];
+ if (not has_boundary_condition[node_id]) {
+ Ar[node_id] = txt * Ar[node_id] * txt + (I - txt);
+
+ has_boundary_condition[node_id] = true;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ // symmetry bc
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
@@ -289,23 +336,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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) {
- } else {
+ if (not has_boundary_condition[node_id]) {
Ar[node_id] = P * Ar[node_id] * P + nxn;
+
+ has_boundary_condition[node_id] = true;
}
}
- } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
-
- parallel_for(
- node_list.size(), PUGS_LAMBDA(size_t i_node) {
- NodeId node_id = node_list[i_node];
-
- Ar[node_id] = identity;
- });
- } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- } else {
- throw UnexpectedError("boundary condition not handled 1");
}
},
boundary_condition);
@@ -329,15 +365,15 @@ 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>
@@ -425,34 +461,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
- NodeValue<int> boundary_node_s(m_mesh.connectivity());
- boundary_node_s.fill(-1);
-
- NodeValue<int> boundary_node_v(m_mesh.connectivity());
- boundary_node_v.fill(-1);
-
+ // pressure bc
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<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId& node_id = node_list[i_node];
- boundary_node_s[node_id] = 1;
- }
- }
- },
- boundary_condition);
- }
-
- 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<VelocityBoundaryCondition, T>) {
- // no changes for velocity conditions
- } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+ 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();
@@ -485,7 +499,91 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// }
// }
// }
- } else if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ }
+ },
+ boundary_condition);
+ }
+
+ NodeValue<bool> has_boundary_condition{m_mesh.connectivity()};
+ has_boundary_condition.fill(false);
+
+ // velocity bc
+ 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<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ if (not has_boundary_condition[node_id]) {
+ has_boundary_condition[node_id] = true;
+ }
+ });
+ }
+ },
+ boundary_condition);
+ }
+
+ if constexpr (Dimension > 1) {
+ // axis bc
+ 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<AxisBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+
+ const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ const Rd& t = bc.direction();
+ const Rdxd txt = tensorProduct(t, t);
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId& node_id = node_list[i_node];
+ if (not has_boundary_condition[node_id]) {
+ const Rdxd inverse_Ar_times_txt = m_inv_Ar[node_id] * txt;
+ 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) * inverse_Ar_times_txt * 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];
+ }
+ }
+ }
+ }
+ }
+ has_boundary_condition[node_id] = true;
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ // symmetry bc
+ 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<SymmetryBoundaryCondition, T>) {
const auto& node_list = bc.nodeList();
const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
@@ -494,13 +592,11 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const Rdxd I = identity;
const Rdxd nxn = tensorProduct(n, n);
const Rdxd P = I - nxn;
- NodeValue<Rdxd> inverse_ArxP{m_mesh.connectivity()};
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) {
- } else {
- inverse_ArxP[node_id] = m_inv_Ar[node_id] * P;
+ if (not has_boundary_condition[node_id]) {
+ const Rdxd inverse_ArxP = m_inv_Ar[node_id] * P;
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);
@@ -517,7 +613,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] * m_Djr(id_cell_j, node_nb_in_j);
+ m_Ajr(id_cell_i, node_nb_in_i) * inverse_ArxP * 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];
}
@@ -525,10 +621,9 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
}
+ has_boundary_condition[node_id] = true;
}
}
- } else {
- throw UnexpectedError("boundary condition not handled 2");
}
},
boundary_condition);
@@ -539,7 +634,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
Vector<double>
- _getU(const DiscreteScalarFunction& p, const DiscreteVectorFunction& u)
+ _getU(const CellValue<const double>& p, const CellValue<const Rd>& u)
{
Vector<double> Un{(Dimension + 1) * m_number_of_implicit_cells};
@@ -565,350 +660,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
Vector<double>
- _getGradJ(const Vector<double>& Un,
- const Vector<double>& Uk,
- const DiscreteScalarFunction& p,
- const DiscreteScalarFunction& pi,
- const DiscreteVectorFunction& u,
- const double dt)
- {
- count_getF++;
- static bool getF_is_started = false;
- if (not getF_is_started) {
- getF_is_started = true;
- getF_t.stop();
- }
-
- getF_t.start();
- 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;
- }();
-
- getF_t.pause();
-
- Vector<double> computed_gradJ{(Dimension + 1) * m_number_of_implicit_cells};
- 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();
-
- 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_Djrpj = zero;
- Rd sum_Ajruj = zero;
- Rd Uk_i = zero;
- Rd Uk_j = zero;
- 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 size_t i_index_p = mapP(id_cell_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);
- Uk_i[i_dimension] = Uk[i_index_u];
- }
- sum_Ajruj += m_Ajr(id_cell_i, node_nb_in_i) * Uk_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 size_t j_index_p = mapP(id_cell_j);
-
- 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);
- Uk_j[i_dimension] = Uk[j_index_u];
- }
- 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_Ajruj +
- m_Ajr(id_cell_j, node_nb_in_j) * Uk_j)[i_dimension];
- }
- }
- }
- }
- }
-
- const double inv_dt = 1. / dt;
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- if (m_is_implicit_cell[cell_id]) {
- size_t j_index_p = mapP(cell_id);
- computed_gradJ[j_index_p] += (inv_dt * m_Mj[cell_id]) * (m_tau_iter[cell_id] - m_tau[cell_id]);
- for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
- size_t j_index_u = mapU(i_dimension, cell_id);
- computed_gradJ[j_index_u] += (inv_dt * m_Mj[cell_id]) * (Uk[j_index_u] - Un[j_index_u]);
- }
- }
- }
-
- for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
- const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
- 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_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_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]];
- }
- }
- for (size_t j_cell = 0; j_cell < node_cells.size(); ++j_cell) {
- const CellId id_cell_j = node_cells[j_cell];
- 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(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_Djrpj_exp +
- m_Ajr(id_cell_j, node_nb_in_j) * m_inv_Ar[node_id] * sum_Ajruj_exp)[i_dimension];
- }
- }
- }
- }
- }
- }
-
- NodeValue<int> boundary_node_s(m_mesh.connectivity());
- boundary_node_s.fill(-1);
-
- NodeValue<int> boundary_node_v(m_mesh.connectivity());
- boundary_node_v.fill(-1);
-
- 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<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId& node_id = node_list[i_node];
- boundary_node_s[node_id] = 1;
- }
- }
- },
- boundary_condition);
- }
-
- for (const auto& boundary_condition : m_boundary_condition_list) {
- std::visit(
- [&](auto&& bc) {
- using T = std::decay_t<decltype(bc)>;
-
- // else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- // const auto& node_list = bc.faceList();
- // const auto& value_list = bc.valueList();
-
- // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(m_mesh);
- // const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
- // const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
- // 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];
- // // 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];
-
- // 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) {
- // const CellId id_cell_j = node_cell[cell_j];
- // if (m_is_implicit_cell[id_cell_j]) {
- // 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(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_Djrpext)[i_dimension];
- // }
- // }
- // }
- // }
- // }
- 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();
- const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
- const Rd& n = bc.outgoingNormal();
- const Rdxd I = identity;
- const Rdxd nxn = tensorProduct(n, n);
- const Rdxd P = I - nxn;
- NodeValue<Rdxd> inverse_ArxP{m_mesh.connectivity()};
- Rd Uk_j = zero;
-
- 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) {
- } else {
- const auto& node_cell = node_to_cell_matrix[node_id];
-
- inverse_ArxP[node_id] = m_inv_Ar[node_id] * P;
- Rd sum_Djrpj = zero;
- Rd sum_Ajruj = 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);
-
- 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 size_t i_index_p = mapP(id_cell_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);
- Uk_i[i_dimension] = Uk[i_index_u];
- }
- sum_Ajruj += m_Ajr(id_cell_i, node_nb_in_i) * Uk_i;
- }
- }
-
- 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_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]];
- }
- }
-
- 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];
-
- int index_p = mapP(id_cell_j);
- computed_gradJ[index_p] +=
- 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);
- Uk_j[i_dimension] = Uk[j_index_u];
- }
-
- 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) * 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_Djrpj_exp +
- m_Ajr(id_cell_j, node_nb_in_j) * inverse_ArxP[node_id] * sum_Ajruj_exp)[i_dimension];
- }
- }
- }
- }
- }
- } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- const auto& value_list = bc.valueList();
-
- 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 (boundary_node_v[node_id] != -1) {
- } else {
- boundary_node_v[node_id] = 1;
- // if (m_is_implicit_node[node_id]) {
- const auto& node_cell = node_to_cell_matrix[node_id];
- Rd Uk_i = zero;
- 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]) {
- 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];
- }
- }
- }
-
- 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 i_index_p = mapP(id_cell_i);
- 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);
- computed_gradJ[i_index_u] +=
- (m_Ajr(id_cell_i, node_nb_in_i) * (Uk_i - value_list[i_node]))[i_dimension];
- }
- }
- }
- }
- }
- } else if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
- } else {
- throw UnexpectedError("boundary condition not handled 3");
- }
- },
- boundary_condition);
- }
- return computed_gradJ;
- }
-
- Vector<double>
- _getF(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)
- {
- Vector<double> gradJ = this->_getGradJ(Un, Uk, p, pi, u, dt);
-
- Vector<double> AU = A * Uk;
-
- Vector<double> F = gradJ - AU;
-
- 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)
+ _getF(const Vector<double>& Un, const Vector<double>& Uk, const DiscreteScalarFunction& pi, const double dt)
{
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -1038,6 +790,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const int i_index_p = mapP(cell_id);
Hess_J(i_index_p, i_index_p) =
(inv_dt * m_Mj[cell_id]) * m_inv_gamma[cell_id] * m_tau_iter[cell_id] / (-Uk[i_index_p] + pi[cell_id]);
+
for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
const int i_index_u = mapU(i_dimension, cell_id);
Hess_J(i_index_u, i_index_u) = inv_dt * m_Mj[cell_id];
@@ -1047,12 +800,29 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
getGradF_t.start();
for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++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);
+
+ 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(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, cell_id_j);
+ for (size_t j_dimension = 0; j_dimension < Dimension; ++j_dimension) {
+ const int j_index_u = mapU(j_dimension, cell_id_j);
+ Hess_J(i_index_u, j_index_u) += Ajr(i_dimension, j_dimension);
+ }
+ }
+ }
+ }
+
if (not is_boundary_node[node_id]) {
const auto& inv_Ar = m_inv_Ar[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);
-
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]) {
@@ -1066,24 +836,11 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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(invArDjr, m_Djr(cell_id_i, node_nb_in_cell_i));
+ Hess_J(i_index_p, j_index_p) += dot(m_Djr(cell_id_i, node_nb_in_cell_i), invArDjr);
}
}
}
}
- 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(cell_id_j, node_nb_in_j_cell);
-
- 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 j_cell = 0; j_cell < node_cell.size(); ++j_cell) {
const CellId cell_id_j = node_cell[j_cell];
@@ -1102,43 +859,186 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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 index_u_i = mapU(j_dimension, cell_id_i);
- Hess_J(index_u_j, index_u_i) -= Ajr_invAr_Air(j_dimension, j_dimension);
+ for (size_t i_dimension = 0; i_dimension < Dimension; ++i_dimension) {
+ const int index_u_i = mapU(i_dimension, cell_id_i);
+ Hess_J(index_u_j, index_u_i) -= Ajr_invAr_Air(j_dimension, i_dimension);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ getGradF_t.pause();
+
+ // presure bc
+ 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>) {
+ // 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(m_Djr(id_cell_i, node_nb_in_i), m_inv_Ar[node_id] * m_Djr(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);
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ }
+ },
+ boundary_condition);
+ }
+
+ NodeValue<bool> has_boundary_condition{m_mesh.connectivity()};
+ has_boundary_condition.fill(false);
+
+ // velocity bc
+ 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<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ has_boundary_condition[node_id] = true;
+ }
+ }
+ },
+ boundary_condition);
+ }
+
+ if constexpr (Dimension > 1) {
+ // axis bc
+ 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<AxisBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+
+ const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ const Rd& t = bc.direction();
+ const Rdxd txt = tensorProduct(t, t);
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId& node_id = node_list[i_node];
+
+ if (not has_boundary_condition[node_id]) {
+ const Rdxd inverse_Ar_times_txt = m_inv_Ar[node_id] * txt;
+
+ 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];
+
+ 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 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) * inverse_Ar_times_txt *
+ m_Ajr(id_cell_j, node_nb_in_j))(i_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 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(m_Djr(id_cell_i, node_nb_in_i), inverse_Ar_times_txt * m_Djr(id_cell_j, node_nb_in_j));
+ }
+ }
+ }
}
+ has_boundary_condition[node_id] = true;
}
}
}
- }
- }
+ },
+ boundary_condition);
}
}
- getGradF_t.pause();
-
- // 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);
-
- NodeValue<int> boundary_node_v(m_mesh.connectivity());
- boundary_node_v.fill(-1);
-
- 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<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
- const NodeId& node_id = node_list[i_node];
- boundary_node_s[node_id] = 1;
- }
- }
- },
- boundary_condition);
- }
+ // symmetry bc
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
@@ -1152,30 +1052,15 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
const Rdxd I = identity;
const Rdxd nxn = tensorProduct(n, n);
const Rdxd P = I - nxn;
- NodeValue<Rdxd> inverse_ArxP{m_mesh.connectivity()};
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) {
- } else {
- inverse_ArxP[node_id] = m_inv_Ar[node_id] * P;
+ if (not has_boundary_condition[node_id]) {
+ const Rdxd inverse_ArxP = m_inv_Ar[node_id] * P;
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];
-
- 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_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]) {
@@ -1192,7 +1077,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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) * inverse_ArxP[node_id] *
+ Hess_J(i_index_u, j_index_u) += (-m_Ajr(id_cell_i, node_nb_in_i) * inverse_ArxP *
m_Ajr(id_cell_j, node_nb_in_j))(i_dimension, j_dimension);
}
}
@@ -1215,103 +1100,14 @@ 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(m_Djr(id_cell_i, node_nb_in_i), inverse_ArxP[node_id] * m_Djr(id_cell_j, node_nb_in_j));
+ dot(m_Djr(id_cell_i, node_nb_in_i), inverse_ArxP * m_Djr(id_cell_j, node_nb_in_j));
}
}
}
}
+ has_boundary_condition[node_id] = true;
}
}
- } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- 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 (boundary_node_v[node_id] != -1) {
- } else {
- boundary_node_v[node_id] = 1;
-
- // if (m_is_implicit_node[node_id]) {
- const auto& node_cell = node_to_cell_matrix[node_id];
- 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 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) {
- int i_index_u = mapU(i_dimension, id_cell_j);
-
- Hess_J(i_index_u, i_index_u) +=
- m_Ajr(id_cell_j, node_local_number_in_its_cells[cell_j])(i_dimension, i_dimension);
- }
- }
- }
- }
- }
- } 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(m_Djr(id_cell_i, node_nb_in_i), m_inv_Ar[node_id] * m_Djr(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");
}
},
boundary_condition);
@@ -1337,7 +1133,25 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CRSMatrix Hess_J_crs = Hess_J.getCRSMatrix();
+ // std::cout << "Hess J = " << Hess_J_crs << '\n';
+
CRSMatrix gradient_f = Hess_J_crs - A;
+
+ // Vector<double> x(gradient_f.numberOfRows());
+ // Vector<double> y(gradient_f.numberOfRows());
+ // double max_err = 0;
+ // for (size_t i = 0; i < gradient_f.numberOfRows(); ++i) {
+ // x.fill(0);
+ // x[i] = 1;
+ // for (size_t j = 0; j < gradient_f.numberOfRows(); ++j) {
+ // y.fill(0);
+ // y[j] = 1;
+
+ // max_err = std::max(max_err, std::abs(dot(x, gradient_f * y) - dot(y, gradient_f * x)));
+ // }
+ // }
+ // std::cout << "MAX NON SYM=" << rang::fgB::cyan << max_err << rang::fg::reset << '\n';
+
HJ_A_t.pause();
return gradient_f;
@@ -1353,6 +1167,18 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
bool is_valid_boundary_condition = true;
switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::axis: {
+ if constexpr (Dimension == 1) {
+ throw NormalError("Axis boundary condition makes no sense in dimension 1");
+ } else {
+ const AxisBoundaryConditionDescriptor& axis_bc_descriptor =
+ dynamic_cast<const AxisBoundaryConditionDescriptor&>(*bc_descriptor);
+
+ bc_list.push_back(
+ AxisBoundaryCondition{getMeshLineNodeBoundary(*mesh, axis_bc_descriptor.boundaryDescriptor())});
+ }
+ break;
+ }
case IBoundaryConditionDescriptor::Type::symmetry: {
const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -1480,12 +1306,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
return implicit_node_index;
}();
- std::cout << "building A: ";
+ std::cout << "building A: " << std::flush;
const CRSMatrix A = this->_getA().getCRSMatrix();
- std::cout << "done\n";
+ std::cout << "done\n" << std::flush;
- Vector<double> Un = this->_getU(p, u);
- Vector<double> Uk = copy(Un);
+ Vector<double> Un = this->_getU(p.cellValues(), u.cellValues());
+ Vector<double> Uk = this->_getU(m_predicted_p, m_predicted_u);
int nb_iter = 0;
double norm_inf_sol;
@@ -1501,17 +1327,12 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
if (m_number_of_implicit_cells > 0) {
do {
nb_iter++;
- // 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);
- std::cout << " done\n";
+ Vector<double> f = this->_getF(Un, Uk, pi, dt);
- // 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';
+ std::cout << "building gradf: " << std::flush;
+ CRSMatrix<double> gradient_f = this->_getGradientF(A, Uk, pi, dt);
+ std::cout << "done\n" << std::flush;
auto l2Norm = [](const Vector<double>& x) {
double sum2 = 0;
@@ -1532,15 +1353,24 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
solver_t.start();
LinearSolver solver;
+ std::cout << "solving linear system: " << std::flush;
solver.solveLocalSystem(gradient_f, sol, f);
+ std::cout << "done\n" << std::flush;
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 - 0.2 * sol;
+ double theta = 1;
+
+ for (CellId cell_id = 0; cell_id < m_number_of_implicit_cells; ++cell_id) {
+ size_t k = mapP(cell_id);
+ if (-Uk[k] + theta * sol[k] < -pi[cell_id]) {
+ double new_theta = 0.5 * (Uk[k] - pi[cell_id]) / sol[k];
+ std::cout << "theta: " << theta << " -> " << new_theta << '\n';
+ theta = new_theta;
+ }
+ }
+
+ Vector<double> U_next = Uk - theta * sol;
Array<const double> abs_sol = [&]() {
Array<double> compute_abs_sol{sol.size()};
@@ -1565,21 +1395,21 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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) {
- // std::cout << "HAS NEG PRESSURE, RELAXING\n";
- // Uk -= 0.01 * sol;
- // } else {
- // Uk -= sol;
- // }
Uk = U_next;
+ // if ((count_newton == 0) and (count_Djr == 0) and (count_timestep == 0)) {
+ // auto newt_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), m_mesh.xr());
+ // double pseudo_time = 1E-8 * count_newton + 1E-3 * count_Djr + count_timestep;
+
+ // std::vector<std::shared_ptr<const INamedDiscreteData>> output_list;
+ // output_list.push_back(
+ // std::make_shared<NamedItemValueVariant>(std::make_shared<ItemValueVariant>(m_predicted_p),
+ // "predicted_p"));
+
+ // vtk_writer.writeOnMeshForced(newt_mesh, output_list, pseudo_time);
+ // }
+
m_predicted_u = [&] {
CellValue<Rd> predicted_u = copy(u.cellValues());
for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
@@ -1590,7 +1420,6 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
predicted_u[cell_id] = vector_u;
}
- // std::cout << "u[" << cell_id << "]=" << predicted_u[cell_id] << '\n';
}
return predicted_u;
}();
@@ -1601,12 +1430,32 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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));
+ NodeValue<const Rd> ur = _getUr();
+
+ NodeValue<Rd> newton_xr{m_mesh.connectivity()};
+ for (NodeId node_id = 0; node_id < m_mesh.numberOfNodes(); ++node_id) {
+ newton_xr[node_id] = m_mesh.xr()[node_id] + dt * ur[node_id];
+ }
+
+ auto newt_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), newton_xr);
+ ++count_newton;
+ // double pseudo_time = 1E-4 * count_newton + 1E-2 * count_Djr + count_timestep;
+
+ // std::vector<std::shared_ptr<const INamedDiscreteData>> output_list;
+ // output_list.push_back(
+ // std::make_shared<NamedItemValueVariant>(std::make_shared<ItemValueVariant>(m_predicted_p), "predicted_p"));
+
+ // vtk_writer.writeOnMeshForced(newt_mesh, output_list, pseudo_time);
+
+ 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';
+ }
+ } while ((norm_inf_sol > 1e-12 * norm_inf_Un) and (nb_iter < 1000));
}
for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
@@ -1640,6 +1489,32 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
b[r] = br;
});
+ NodeValue<bool> has_boundary_condition{m_mesh.connectivity()};
+ has_boundary_condition.fill(false);
+
+ // velocity bc
+ 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<VelocityBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+ const auto& value_list = bc.valueList();
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& value = value_list[i_node];
+ b[node_id] = value;
+
+ has_boundary_condition[node_id] = true;
+ });
+ }
+ },
+ boundary_condition);
+ }
+
+ // pressure bc
for (const auto& boundary_condition : m_boundary_condition_list) {
std::visit(
[&](auto&& bc) {
@@ -1660,8 +1535,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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);
+ const CellId node_cell_id = node_cell_list[0];
+ const 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);
});
@@ -1680,8 +1555,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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 CellId face_cell_id = face_cell_list[0];
+ const 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;
@@ -1693,33 +1568,57 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
}
}
}
- } 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];
- });
+ }
+ },
+ boundary_condition);
+ }
- } else if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
- const auto& node_list = bc.nodeList();
- const auto& value_list = bc.valueList();
+ if constexpr (Dimension > 1) {
+ // axis bc
+ 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<AxisBoundaryCondition, T>) {
+ const Rd& t = bc.direction();
+ const Rdxd txt = tensorProduct(t, t);
+
+ const auto& node_list = bc.nodeList();
+ parallel_for(
+ bc.numberOfNodes(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ if (not has_boundary_condition[node_id]) {
+ b[node_id] = txt * b[node_id];
+
+ has_boundary_condition[node_id] = true;
+ }
+ });
+ }
+ },
+ boundary_condition);
+ }
+ }
+ // symmetry bc
+ 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<SymmetryBoundaryCondition, T>) {
+ const Rd& n = bc.outgoingNormal();
+ const Rdxd I = identity;
+ const Rdxd nxn = tensorProduct(n, n);
+ const Rdxd P = I - nxn;
+ const auto& node_list = bc.nodeList();
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;
+ bc.numberOfNodes(), PUGS_LAMBDA(const size_t i_node) {
+ const NodeId node_id = node_list[i_node];
+ if (not has_boundary_condition[node_id]) {
+ b[node_id] = P * b[node_id];
+
+ has_boundary_condition[node_id] = true;
+ }
});
- } else {
- throw UnexpectedError("boundary condition not handled 5");
}
},
boundary_condition);
@@ -1835,6 +1734,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
const NodeValuePerCell<const Rd> Cjr_n = mesh_data.Cjr();
+ count_Djr = 0;
+
m_Djr = copy(Cjr_n);
CellValue<double> new_rho = copy(rho.cellValues());
@@ -1892,9 +1793,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
int number_iter = 0;
do {
number_iter++;
- new_rho = copy(rho.cellValues());
- new_u = copy(u.cellValues());
- new_E = copy(E.cellValues());
+
+ count_newton = 0;
this->_computeGradJU_AU(u, p, pi, dt);
@@ -1912,6 +1812,10 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+ new_rho = copy(rho.cellValues());
+ new_u = copy(u.cellValues());
+ new_E = copy(E.cellValues());
+
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_nodes = cell_to_node_matrix[j];
@@ -1930,19 +1834,33 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// 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_Vj[j] = 0;
- auto cell_nodes = cell_to_node_matrix[j];
+ 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];
- 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;
- });
+ 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;
+ });
+
+ double m = min(new_Vj);
+ if (m < 0) {
+ std::cout << "negative volume\n";
+ 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);
parallel_for(
mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
@@ -2011,10 +1929,7 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
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])),
+ 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;
@@ -2077,12 +1992,31 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
});
m_Djr = Djr;
+
+ 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));
+ }
+ }
+
+ std::cout << rang::fgB::yellow << "Max volume err = " << max_err << rang::fg::reset << '\n';
}
if constexpr (Dimension > 1) {
- std::cout << "number_iter_Djr=" << number_iter << " max_tau_error=" << max_tau_error << '\n';
+ std::cout << rang::fgB::magenta << "number_iter_Djr=" << number_iter << " max_tau_error=" << max_tau_error
+ << rang::fg::reset << '\n';
}
+ ++count_Djr;
+
// 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';
@@ -2093,6 +2027,8 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver final
// << " HJ_A=" << HJ_A_t << '\n';
// std::cout << "solver=" << solver_t << " total= " << implicit_t << '\n';
+ ++count_timestep;
+
return {new_mesh, std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction{new_mesh, new_rho}),
std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction{new_mesh, new_u}),
std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction{new_mesh, new_E})};
@@ -2308,6 +2244,51 @@ class ImplicitAcousticSolverHandler::ImplicitAcousticSolver<Dimension>::Symmetry
~SymmetryBoundaryCondition() = default;
};
+template <size_t Dimension>
+class ImplicitAcousticSolverHandler::ImplicitAcousticSolver<Dimension>::AxisBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
+
+ public:
+ const Rd&
+ direction() const
+ {
+ return m_mesh_line_node_boundary.direction();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_line_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_line_node_boundary.nodeList();
+ }
+
+ AxisBoundaryCondition(const MeshLineNodeBoundary<Dimension>& mesh_line_node_boundary)
+ : m_mesh_line_node_boundary(mesh_line_node_boundary)
+ {
+ ;
+ }
+
+ ~AxisBoundaryCondition() = default;
+};
+
+template <>
+class ImplicitAcousticSolverHandler::ImplicitAcousticSolver<1>::AxisBoundaryCondition
+{
+ public:
+ AxisBoundaryCondition() = default;
+ ~AxisBoundaryCondition() = default;
+};
+
ImplicitAcousticSolverHandler::ImplicitAcousticSolverHandler(
const SolverType solver_type,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,