LocalDtAcousticSolver.cpp

#include <scheme/LocalDtAcousticSolver.hpp>

#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/ItemValueVariant.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
#include <scheme/FixedBoundaryConditionDescriptor.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>

#include <variant>
#include <vector>

template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
  : public LocalDtAcousticSolverHandler::ILocalDtAcousticSolver
{
 private:
  using Rdxd = TinyMatrix<Dimension>;
  using Rd   = TinyVector<Dimension>;

  using MeshType     = Mesh<Connectivity<Dimension>>;
  using MeshDataType = MeshData<Dimension>;

  using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, const double>;
  using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, const Rd>;

  class FixedBoundaryCondition;
  class PressureBoundaryCondition;
  class SymmetryBoundaryCondition;
  class VelocityBoundaryCondition;
  class ExternalFSIVelocityBoundaryCondition;

  using BoundaryCondition = std::variant<FixedBoundaryCondition,
                                         PressureBoundaryCondition,
                                         SymmetryBoundaryCondition,
                                         VelocityBoundaryCondition,
                                         ExternalFSIVelocityBoundaryCondition>;

  using BoundaryConditionList = std::vector<BoundaryCondition>;

  NodeValuePerCell<const Rdxd>
  _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    const NodeValuePerCell<const Rd> njr = mesh_data.njr();

    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const size_t& nb_nodes = Ajr.numberOfSubValues(j);
        const double& rhoc_j   = rhoc[j];
        for (size_t r = 0; r < nb_nodes; ++r) {
          Ajr(j, r) = tensorProduct(rhoc_j * Cjr(j, r), njr(j, r));
        }
      });

    return Ajr;
  }

  NodeValuePerCell<const Rdxd>
  _computeEucclhydAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
    const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();

    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();

    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};

    parallel_for(
      Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });

    parallel_for(
      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 double& rho_c = rhoc[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();
          };

          for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
            const size_t R = local_node_number_in_cell(face_nodes[rl]);
            Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
          }
        }
      });

    return Ajr;
  }

  NodeValuePerCell<const Rdxd>
  _computeAjr(const SolverType& solver_type, const MeshType& mesh, const DiscreteScalarFunction& rhoc) const
  {
    if constexpr (Dimension == 1) {
      return _computeGlaceAjr(mesh, rhoc);
    } else {
      switch (solver_type) {
      case SolverType::Glace: {
        return _computeGlaceAjr(mesh, rhoc);
      }
      case SolverType::Eucclhyd: {
        return _computeEucclhydAjr(mesh, rhoc);
      }
      default: {
        throw UnexpectedError("invalid solver type");
      }
      }
    }
  }

  NodeValue<Rdxd>
  _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
  {
    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValue<Rdxd> Ar{mesh.connectivity()};

    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        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);

        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);
        }
        Ar[r] = sum;
      });

    return Ar;
  }

  NodeValue<Rd>
  _computeBr(const Mesh<Connectivity<Dimension>>& mesh,
             const NodeValuePerCell<const Rdxd>& Ajr,
             const DiscreteVectorFunction& u,
             const DiscreteScalarFunction& p) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();

    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

    NodeValue<Rd> b{mesh.connectivity()};

    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
        const auto& node_to_cell                   = node_to_cell_matrix[r];
        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);

        Rd br = zero;
        for (size_t j = 0; j < node_to_cell.size(); ++j) {
          const CellId J       = node_to_cell[j];
          const unsigned int R = node_local_number_in_its_cells[j];
          br += Ajr(J, R) * u[J] + p[J] * Cjr(J, R);
        }

        b[r] = br;
      });

    return b;
  }

  BoundaryConditionList
  _getBCList(const MeshType& mesh,
             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
  {
    BoundaryConditionList bc_list;

    for (const auto& bc_descriptor : bc_descriptor_list) {
      bool is_valid_boundary_condition = true;

      switch (bc_descriptor->type()) {
      case IBoundaryConditionDescriptor::Type::symmetry: {
        bc_list.emplace_back(
          SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        break;
      }
      case IBoundaryConditionDescriptor::Type::fixed: {
        bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
        break;
      }
      case IBoundaryConditionDescriptor::Type::external: {
        if constexpr (Dimension == 1) {
          const ExternalBoundaryConditionDescriptor& extern_bc_descriptor =
            dynamic_cast<const ExternalBoundaryConditionDescriptor&>(*bc_descriptor);
          bc_list.emplace_back(
            ExternalFSIVelocityBoundaryCondition(mesh, getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
                                                 extern_bc_descriptor.socket()));
        } else {
          throw NotImplementedError("external FSI velocity not implemented for dimension > 1");
        }
        break;
      }
      case IBoundaryConditionDescriptor::Type::dirichlet: {
        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
        if (dirichlet_bc_descriptor.name() == "velocity") {
          MeshNodeBoundary<Dimension> mesh_node_boundary =
            getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());

          Array<const Rd> value_list =
            InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
                                                                               mesh.xr(),
                                                                               mesh_node_boundary.nodeList());

          bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
        } else if (dirichlet_bc_descriptor.name() == "pressure") {
          const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();

          if constexpr (Dimension == 1) {
            MeshNodeBoundary<Dimension> mesh_node_boundary =
              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());

            Array<const double> node_values =
              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
                                                                                     mesh_node_boundary.nodeList());

            bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
          } else {
            MeshFaceBoundary<Dimension> mesh_face_boundary =
              getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());

            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
            Array<const double> face_values =
              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
                                                                                     mesh_face_boundary.faceList());
            bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
          }

        } else {
          is_valid_boundary_condition = false;
        }
        break;
      }
      default: {
        is_valid_boundary_condition = false;
      }
      }
      if (not is_valid_boundary_condition) {
        std::ostringstream error_msg;
        error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
        throw NormalError(error_msg.str());
      }
    }

    return bc_list;
  }

  void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
  void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
  void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
  void _applyExternalVelocityBC(const BoundaryConditionList& bc_list,
                                const DiscreteScalarFunction& p,
                                NodeValue<Rdxd>& Ar,
                                NodeValue<Rd>& br) const;

  void
  _applyBoundaryConditions(const BoundaryConditionList& bc_list,
                           const MeshType& mesh,
                           NodeValue<Rdxd>& Ar,
                           NodeValue<Rd>& br) const
  {
    this->_applyPressureBC(bc_list, mesh, br);
    this->_applySymmetryBC(bc_list, Ar, br);
    this->_applyVelocityBC(bc_list, Ar, br);
  }

  NodeValue<const Rd>
  _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
  {
    NodeValue<Rd> u{mesh.connectivity()};
    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });

    return u;
  }

  NodeValuePerCell<Rd>
  _computeFjr(const MeshType& mesh,
              const NodeValuePerCell<const Rdxd>& Ajr,
              const NodeValue<const Rd>& ur,
              const DiscreteVectorFunction& u,
              const DiscreteScalarFunction& p) const
  {
    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);

    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();

    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    NodeValuePerCell<Rd> F{mesh.connectivity()};
    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_nodes = cell_to_node_matrix[j];

        for (size_t r = 0; r < cell_nodes.size(); ++r) {
          F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r);
        }
      });

    return F;
  }

 public:
  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
  compute_fluxes(const SolverType& solver_type,
                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
                 const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
  {
    std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
    if (not i_mesh) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
      throw NormalError("acoustic solver expects P0 functions");
    }

    const MeshType& mesh              = dynamic_cast<const MeshType&>(*i_mesh);
    const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
    const DiscreteScalarFunction& c   = c_v->get<DiscreteScalarFunction>();
    const DiscreteVectorFunction& u   = u_v->get<DiscreteVectorFunction>();
    const DiscreteScalarFunction& p   = p_v->get<DiscreteScalarFunction>();

    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);

    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);

    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
    this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
    this->_applyExternalVelocityBC(bc_list, p, Ar, br);

    synchronize(Ar);
    synchronize(br);

    NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);

    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                           std::make_shared<const SubItemValuePerItemVariant>(Fjr));
  }

  std::tuple<std::shared_ptr<const IMesh>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>>
  apply_fluxes(const double& dt,
               const MeshType& mesh,
               const DiscreteScalarFunction& rho,
               const DiscreteVectorFunction& u,
               const DiscreteScalarFunction& E,
               const NodeValue<const Rd>& ur,
               const NodeValuePerCell<const Rd>& Fjr) const
  {
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
        (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
      throw NormalError("fluxes are not defined on the same connectivity than the mesh");
    }

    NodeValue<Rd> new_xr = copy(mesh.xr());
    parallel_for(
      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });

    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);

    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();

    CellValue<double> new_rho = copy(rho.cellValues());
    CellValue<Rd> new_u       = copy(u.cellValues());
    CellValue<double> new_E   = copy(E.cellValues());

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
        const auto& cell_nodes = cell_to_node_matrix[j];

        Rd momentum_fluxes   = zero;
        double energy_fluxes = 0;
        for (size_t R = 0; R < cell_nodes.size(); ++R) {
          const NodeId r = cell_nodes[R];
          momentum_fluxes += Fjr(j, R);
          energy_fluxes += dot(Fjr(j, R), ur[r]);
        }
        const double dt_over_Mj = dt / (rho[j] * Vj[j]);
        new_u[j] -= dt_over_Mj * momentum_fluxes;
        new_E[j] -= dt_over_Mj * energy_fluxes;
      });

    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();

    parallel_for(
      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });

    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))};
  }

  std::tuple<std::shared_ptr<const IMesh>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>>
  apply_fluxes(const double& dt,
               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
               const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
               const std::shared_ptr<const ItemValueVariant>& ur,
               const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
  {
    std::shared_ptr i_mesh = getCommonMesh({rho_v, u_v, E_v});
    if (not i_mesh) {
      throw NormalError("discrete functions are not defined on the same mesh");
    }

    if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
      throw NormalError("acoustic solver expects P0 functions");
    }

    return this->apply_fluxes(dt,                                       //
                              dynamic_cast<const MeshType&>(*i_mesh),   //
                              rho_v->get<DiscreteScalarFunction>(),     //
                              u_v->get<DiscreteVectorFunction>(),       //
                              E_v->get<DiscreteScalarFunction>(),       //
                              ur->get<NodeValue<const Rd>>(),           //
                              Fjr->get<NodeValuePerCell<const Rd>>());
  }

  std::tuple<std::shared_ptr<const IMesh>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>,
             std::shared_ptr<const DiscreteFunctionVariant>>
  apply(const SolverType& solver_type,
        const double& dt,
        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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
  {
    auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
    return apply_fluxes(dt, rho, u, E, ur, Fjr);
  }

  LocalDtAcousticSolver()                        = default;
  LocalDtAcousticSolver(LocalDtAcousticSolver&&) = default;
  ~LocalDtAcousticSolver()                       = default;
};

template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
                                                                                 const MeshType& mesh,
                                                                                 NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_list) {
    std::visit(
      [&](auto&& bc) {
        using T = std::decay_t<decltype(bc)>;
        if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
          MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
          if constexpr (Dimension == 1) {
            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();

            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
            const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();

            const auto& node_list  = bc.nodeList();
            const auto& value_list = bc.valueList();
            parallel_for(
              node_list.size(), PUGS_LAMBDA(size_t i_node) {
                const NodeId node_id       = node_list[i_node];
                const auto& node_cell_list = node_to_cell_matrix[node_id];
                Assert(node_cell_list.size() == 1);

                CellId node_cell_id              = node_cell_list[0];
                size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);

                br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
              });
          } else {
            const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();

            const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
            const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();

            const auto& face_list  = bc.faceList();
            const auto& value_list = bc.valueList();
            for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
              const FaceId face_id       = face_list[i_face];
              const auto& face_cell_list = face_to_cell_matrix[face_id];
              Assert(face_cell_list.size() == 1);

              CellId face_cell_id              = face_cell_list[0];
              size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);

              const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;

              const auto& face_nodes = face_to_node_matrix[face_id];

              for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
                NodeId node_id = face_nodes[i_node];
                br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
              }
            }
          }
        }
      },
      boundary_condition);
  }
}

template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
                                                                                 NodeValue<Rdxd>& Ar,
                                                                                 NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_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 Array<const NodeId>& node_list = bc.nodeList();
          parallel_for(
            bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
              const NodeId r = node_list[r_number];

              Ar[r] = P * Ar[r] * P + nxn;
              br[r] = P * br[r];
            });
        }
      },
      boundary_condition);
  }
}

template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
                                                                                 NodeValue<Rdxd>& Ar,
                                                                                 NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_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(size_t i_node) {
              NodeId node_id    = node_list[i_node];
              const auto& value = value_list[i_node];

              Ar[node_id] = identity;
              br[node_id] = value;
            });
        } else if constexpr (std::is_same_v<FixedBoundaryCondition, 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;
              br[node_id] = zero;
            });
        }
      },
      boundary_condition);
  }
}

template <size_t Dimension>
void
LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyExternalVelocityBC(
  const BoundaryConditionList& bc_list,
  const DiscreteScalarFunction& p,
  NodeValue<Rdxd>& Ar,
  NodeValue<Rd>& br) const
{
  for (const auto& boundary_condition : bc_list) {
    std::visit(
      [&](auto&& bc) {
        using T = std::decay_t<decltype(bc)>;
        if constexpr (std::is_same_v<ExternalFSIVelocityBoundaryCondition, T>) {
          const auto& node_list  = bc.nodeList();
          const auto& value_list = bc.valueList(p);

          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];

              Ar[node_id] = identity;
              br[node_id] = value;
            });
        }
      },
      boundary_condition);
  }
}

template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::FixedBoundaryCondition
{
 private:
  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
    : m_mesh_node_boundary{mesh_node_boundary}
  {}

  ~FixedBoundaryCondition() = default;
};

template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::PressureBoundaryCondition
{
 private:
  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
  const Array<const double> m_value_list;

 public:
  const Array<const FaceId>&
  faceList() const
  {
    return m_mesh_face_boundary.faceList();
  }

  const Array<const double>&
  valueList() const
  {
    return m_value_list;
  }

  PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
                            const Array<const double>& value_list)
    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
  {}

  ~PressureBoundaryCondition() = default;
};

template <>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<1>::PressureBoundaryCondition
{
 private:
  const MeshNodeBoundary<1> m_mesh_node_boundary;
  const Array<const double> m_value_list;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const double>&
  valueList() const
  {
    return m_value_list;
  }

  PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
  {}

  ~PressureBoundaryCondition() = default;
};

template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::ExternalFSIVelocityBoundaryCondition
{
 private:
  const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
  const Array<TinyVector<Dimension>> m_value_list;
  const std::shared_ptr<const Socket> m_socket;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  Array<const TinyVector<Dimension>>
  valueList(const DiscreteScalarFunction& p) const
  {
    if (parallel::size() > 1) {
      throw NotImplementedError("Parallelism is not supported");
    }
    if (m_value_list.size() != 1) {
      throw NotImplementedError("Non connected boundaries are not supported");
    }

    const CellId cell_id = m_node_to_cell_matrix[m_mesh_node_boundary.nodeList()[0]][0];

    write(*m_socket, p[cell_id]);
    read(*m_socket, m_value_list[0]);
    return m_value_list;
  }

  ExternalFSIVelocityBoundaryCondition(const Mesh<Connectivity<Dimension>>& mesh,
                                       const MeshNodeBoundary<Dimension>& mesh_node_boundary,
                                       const std::shared_ptr<const Socket>& socket)
    : m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()},
      m_mesh_node_boundary{mesh_node_boundary},
      m_value_list(mesh_node_boundary.nodeList().size()),
      m_socket{socket}
  {
    if constexpr (Dimension > 1) {
      throw NotImplementedError("external FSI velocity bc in dimension > 1");
    }
  }

  ~ExternalFSIVelocityBoundaryCondition() = default;
};

template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::VelocityBoundaryCondition
{
 private:
  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;

  const Array<const TinyVector<Dimension>> m_value_list;

 public:
  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_node_boundary.nodeList();
  }

  const Array<const TinyVector<Dimension>>&
  valueList() const
  {
    return m_value_list;
  }

  VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
                            const Array<const TinyVector<Dimension>>& value_list)
    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
  {}

  ~VelocityBoundaryCondition() = default;
};

template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::SymmetryBoundaryCondition
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_mesh_flat_node_boundary.outgoingNormal();
  }

  size_t
  numberOfNodes() const
  {
    return m_mesh_flat_node_boundary.nodeList().size();
  }

  const Array<const NodeId>&
  nodeList() const
  {
    return m_mesh_flat_node_boundary.nodeList();
  }

  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
  {
    ;
  }

  ~SymmetryBoundaryCondition() = default;
};

LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh)
{
  if (not i_mesh) {
    throw NormalError("discrete functions are not defined on the same mesh");
  }

  switch (i_mesh->dimension()) {
  case 1: {
    m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<1>>();
    break;
  }
  case 2: {
    m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<2>>();
    break;
  }
  case 3: {
    m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<3>>();
    break;
  }
  default: {
    throw UnexpectedError("invalid mesh dimension");
  }
  }
}