add some corrections on Boundary conditions

[ci skip]

src/scheme/CMakeLists.txt

@@ -6,4 +6,5 @@ add_library(
   DiscreteFunctionInterpoler.cpp
   DiscreteFunctionUtils.cpp
   ImplicitAcousticSolver.cpp
+  UpwindExplicitTrafficSolver.cpp
   PerfectGas.cpp)

src/scheme/UpwindExplicitTrafficSolver.cpp

 #include <scheme/UpwindExplicitTrafficSolver.hpp>
 
+#include <language/utils/FunctionSymbolId.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-
 #include <variant>
 #include <vector>
 
@@ -36,6 +37,11 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
   BoundaryConditionList m_boundary_condition_list;
   const MeshType& m_mesh;
 
+  NodeValue<const Rd> m_ur;
+  NodeValue<const Rd> m_upwind_flux;
+  CellValue<const Rd> m_tauj_carre;
+  CellValue<const Rd> m_Mj;
+
   CellValue<const double>
   _getRho(const DiscreteScalarFunction& tau)
   {
@@ -85,12 +91,13 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
     return inv_Mj;
   }
 
-  const NodeValue<const double> _getUpwindFlux = (const DiscreteScalarFunction& tau)[&]()
+  NodeValue<const Rd>
+  _getUpwindFlux(const DiscreteScalarFunction& tau)
   {
     const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
 
     const auto& node_local_numbers_in_their_cells = m_mesh.connectivity().nodeLocalNumbersInTheirCells();
-    NodeValue<double> computed_flux(m_mesh.connectivity());
+    NodeValue<Rd> computed_flux(m_mesh.connectivity());
     parallel_for(
       m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
         const auto& node_cells = node_to_cell_matrix[r];
@@ -128,7 +135,83 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
     return computed_flux;
   }
 
-  NodeValue<const Rd> _computeUr = (const MeshType& mesh, NodeValue<const double> upwind_flux)
+  BoundaryConditionList
+  _getBCList(const std::shared_ptr<const MeshType>& mesh,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+  {
+    BoundaryConditionList bc_list;
+
+    constexpr ItemType FaceType = [] {
+      if constexpr (Dimension > 1) {
+        return ItemType::face;
+      } else {
+        return ItemType::node;
+      }
+    }();
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      bool is_valid_boundary_condition = true;
+
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::dirichlet: {
+        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+        if (dirichlet_bc_descriptor.name() == "velocity") {
+          for (size_t i_ref_face_list = 0;
+               i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+            const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+            const RefId& ref          = ref_face_list.refId();
+            if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
+              MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+
+              const FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+              const auto& node_list = mesh_node_boundary.nodeList();
+
+              Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+                TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list);
+
+              bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+            }
+          }
+        } else {
+          is_valid_boundary_condition = false;
+        }
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::free: {   // <- Free boundary
+        const FreeBoundaryConditionDescriptor& free_bc_descriptor =
+          dynamic_cast<const FreeBoundaryConditionDescriptor&>(*bc_descriptor);
+
+        for (size_t i_ref_face_list = 0;
+             i_ref_face_list < mesh->connectivity().template numberOfRefItemList<ItemType::node>(); ++i_ref_face_list) {
+          const auto& ref_face_list = mesh->connectivity().template refItemList<ItemType::node>(i_ref_face_list);
+          const RefId& ref          = ref_face_list.refId();
+          if (ref == free_bc_descriptor.boundaryDescriptor()) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+            const auto& node_list = mesh_node_boundary.nodeList();
+
+            bc_list.push_back(FreeBoundaryCondition{node_list});
+          }
+        }
+        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;
+  }
+
+  NodeValue<const Rd>
+  _computeUr(const MeshType& mesh, NodeValue<const double> upwind_flux)
   {
     NodeValue<Rd> computed_ur(mesh.connectivity());
     parallel_for(
@@ -140,32 +223,34 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
   UpwindExplicitTrafficSolver(
     const std::shared_ptr<const MeshType>& p_mesh,
     const DiscreteScalarFunction tau,
-    const std::vector<std::shared_ptr<const IBoundaryConditionDescrptior>>& bc_descriptor_list)
-    : m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+    : m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}, m_mesh{*p_mesh}
   {
     const MeshType& mesh = *p_mesh;
 
     CellValue<double> rhoj   = this->_getRho(tau);
     CellValue<double> uj     = this->_getU(tau);
-    CellValue<double> tauj_carre = this->_getTau_carre(tau);
-    CellValue<double> Mj         = this->_getMj(tau);
     CellValue<double> inv_Mj = this->_getInv_Mj(tau);
 
     m_upwind_flux = this->_getUpwindFlux(tau);
     m_ur          = this->_computeUr(mesh, m_upwind_flux);
+    m_tauj_carre  = this->_getTaucarre(tau);
+    m_Mj          = this->_getMj(tau);
+    ;
   }
 
  public:
   std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
   apply(const std::shared_ptr<const MeshType>& mesh, const DiscreteFunctionP0<Dimension, double>& tau) const
   {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
     double c                = 0;
-    c        = 1. / min(tauj_carre);
+    c                       = 1. / min(m_tauj_carre);
 
     double dt = 0;
-    dt        = 0.5 * min(Mj) / c;
+    dt        = 0.5 * min(m_Mj) / c;
 
-    const NodeValuePerCell<const TinyVector<1, double>> Cjr = mesh_data.Cjr();
+    const NodeValuePerCell<Rd> Cjr = mesh_data.Cjr();
 
     const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
 
@@ -177,7 +262,7 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
 
     CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
 
-    const NodeValuePerCell<const Rd>& Cjr = MeshDataManager::instance().getMeshData(*mesh).Cjr();
+    Cjr = MeshDataManager::instance().getMeshData(*mesh).Cjr();
 
     CellValue<double> new_tau = copy(tau.cellValues());
 
@@ -195,9 +280,6 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
 
     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<DiscreteScalarFunction>(new_mesh, new_tau)};
   }
   std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>>
@@ -228,34 +310,34 @@ class UpwindExplicitTrafficSolverHandler ::UpwindExplicitTrafficSolver final
 };
 
 template <size_t Dimension>
-void
-UpwindExplicitTrafficSolverHandler::UpwindExplicitTrafficSolver<Dimension>::_applyVelocityBC(NodeValue<Rdxd>& Ar,
-                                                                                             NodeValue<Rd>& br)
+class UpwindExplicitTrafficSolverHandler::UpwindExplicitTrafficSolver<Dimension>::VelocityBoundaryCondition
 {
-  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(size_t i_node) {
-              NodeId node_id    = node_list[i_node];
-              const auto& value = value_list[i_node];
+ private:
+  const Array<const TinyVector<Dimension>> m_value_list;
+  const Array<const NodeId> m_node_list;
 
-              Ar[node_id] = identity;
-              br[node_id] = value;
-            });
-        }
-      },
-      boundary_condition);
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_node_list;
   }
+
+  const Array<const TinyVector<Dimension>>&
+  valueList() const
+  {
+    return m_value_list;
   }
 
+  VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
+    : m_value_list{value_list}, m_node_list{node_list}
+  {}
+
+  ~VelocityBoundaryCondition() = default;
+};
+
 template <size_t Dimension>
-class UpwindExplicitTrafficHandler::UpwindExplicitTrafficSolver<Dimension>::FreeBoundaryCondition
+class UpwindExplicitTrafficSolverHandler::UpwindExplicitTrafficSolver<Dimension>::FreeBoundaryCondition
 {
  private:
   const Array<const NodeId> m_node_list;
@@ -272,13 +354,7 @@ class UpwindExplicitTrafficHandler::UpwindExplicitTrafficSolver<Dimension>::Free
   ~FreeBoundaryCondition() = default;
 };
 
-template UpwindExplicitTrafficHandler::UpwindExplicitTrafficSolver<1>::UpwindExplicitTrafficAlgorithm(
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const double&);
-
-UpwindExplicitTrafficHandler::UpwindExplicitTrafficSolverHandler(
+UpwindExplicitTrafficSolverHandler::UpwindExplicitTrafficSolverHandler(
   const std::shared_ptr<const IDiscreteFunction>& tau,
   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
 {
@@ -293,15 +369,18 @@ UpwindExplicitTrafficHandler::UpwindExplicitTrafficSolverHandler(
 
   switch (i_mesh->dimension()) {
   case 1: {
-    m_traffic_solver = std::make_unique<UpwindExplicitTrafficSolver<1>>(i_mesh, tau, bc_descriptor_list);
+    m_upwind_explicit_traffic_solver =
+      std::make_unique<UpwindExplicitTrafficSolver<1>>(i_mesh, tau, bc_descriptor_list);
     break;
   }
   case 2: {
-    m_traffic_solver = std::make_unique<UpwindExplicitTrafficSolver<2>>(i_mesh, tau, bc_descriptor_list);
+    m_upwind_explicit_traffic_solver =
+      std::make_unique<UpwindExplicitTrafficSolver<2>>(i_mesh, tau, bc_descriptor_list);
     break;
   }
   case 3: {
-    m_traffic_solver = std::make_unique<UpwindExplicitTrafficSolver<3>>(i_mesh, tau, bc_descriptor_list);
+    m_upwind_explicit_traffic_solver =
+      std::make_unique<UpwindExplicitTrafficSolver<3>>(i_mesh, tau, bc_descriptor_list);
     break;
   }
   default: {

src/scheme/UpwindExplicitTrafficSolver.hpp

@@ -5,7 +5,7 @@
 #include <tuple>
 #include <vector>
 
-class IDiscreteFUnction;
+class IDiscreteFunction;
 class IBoundaryConditionDescriptor;
 class IMesh;