feat[femutils]: start working on applying dirichlet via elimination f…

…or bsr format

femutils/AlephDoFLinearSystem.cc

@@ -21,6 +21,7 @@
 
 #include <arcane/aleph/AlephTypesSolver.h>
 #include <arcane/aleph/Aleph.h>
+#include <arccore/base/NotImplementedException.h>
 
 #include "FemUtils.h"
 #include "IDoFLinearSystemFactory.h"
@@ -245,6 +246,11 @@ class AlephDoFLinearSystemImpl
     info() << "EliminateRowColumn row=" << row.localId() << " v=" << value;
   }
 
+  void setEliminationArrays(VariableDoFByte& dof_elimination_info, VariableDoFReal& dof_elimination_value) override
+  {
+    ARCANE_THROW(NotImplementedException, "");
+  };
+
   void solve() override
   {
     UniqueArray<Real> aleph_result;

femutils/BSRFormat.h

@@ -53,6 +53,7 @@
 #include "arcane/accelerator/VariableViews.h"
 #include <arcane/accelerator/Atomic.h>
 #include <arcane/accelerator/Scan.h>
+#include <type_traits>
 
 #include "DoFLinearSystem.h"
 #include "CsrFormatMatrix.h"
@@ -824,8 +825,8 @@ class BSRFormat : public TraceAccessor
   /*---------------------------------------------------------------------------*/
   /*---------------------------------------------------------------------------*/
 
-  template <bool IS_WEAK, bool USE_CSR_IN_LINEAR_SYSTEM>
-  void applyDirichletByPenalty(Real penalty, NumArray<Real, MDDim1>& rhs_vect, std::array<VariableNodeByte, NB_DOF>& u_dirichlet_var_arr, VariableNodeReal u)
+  template <bool IS_WEAK, bool USE_CSR_IN_LINEAR_SYSTEM, typename DataType>
+  void applyDirichletViaPenalty(Real penalty, NumArray<Real, MDDim1>& rhs_vect, std::array<VariableNodeByte, NB_DOF>& u_dirichlet_var_arr, MeshVariableScalarRefT<Node, DataType> u)
   {
     auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
     auto nb_row = m_bsr_matrix.nbRow();
@@ -876,54 +877,62 @@ class BSRFormat : public TraceAccessor
   /*---------------------------------------------------------------------------*/
 
   template <bool ELIMINATE_BY_ROW_N_COLUMN, typename DataType>
-  void applyDirichletByElimination(std::array<VariableNodeByte, NB_DOF>& u_dirichlet_var_arr, MeshVariableScalarRefT<Node, DataType> u, DoFLinearSystem& linear_system)
+  void applyDirichletViaElimination(std::array<VariableNodeByte, NB_DOF>& u_dirichlet_var_arr, MeshVariableScalarRefT<Node, DataType> u, DoFLinearSystem* linear_system)
   {
     auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
     auto command = makeCommand(m_queue);
     auto in_u = viewIn(command, u);
+    VariableDoFByte m_dof_elimination_info(VariableBuildInfo(m_dofs_on_nodes.dofFamily(), "BSRFormatDoFEliminationInfo"));
+    VariableDoFReal m_dof_elimination_value(VariableBuildInfo(m_dofs_on_nodes.dofFamily(), "BSRFormatDoFEliminationValue"));
+    auto in_out_dof_elimination_info = viewInOut(command, m_dof_elimination_info);
+    auto in_out_dof_elimination_value = viewInOut(command, m_dof_elimination_value);
 
     NumArray<Accelerator::ItemVariableScalarInViewT<Node, Byte>, MDDim1> in_u_dirichlet_view_arr;
     in_u_dirichlet_view_arr.resize(NB_DOF);
     for (auto i = 0; i < NB_DOF; ++i)
       in_u_dirichlet_view_arr[i] = viewIn(command, u_dirichlet_var_arr[i]);
     auto in_u_dirichlet_arr = viewIn(command, in_u_dirichlet_view_arr);
 
-    /*auto eliminateInLinearSystem = [&linear_system] ARCCORE_HOST_DEVICE(DoFLocalId dof_lid, Real u_dirichlet) {
-      if constexpr (ELIMINATE_BY_ROW_N_COLUMN)
-        linear_system.eliminateRowColumn(dof_lid, u_dirichlet);
-      else
-        linear_system.eliminateRow(dof_lid, u_dirichlet);
-    };*/
-
-    command << RUNCOMMAND_ENUMERATE(NodeLocalId, node_lid, m_mesh->ownNodes())
-    {
-      for (auto i = 0; i < NB_DOF; ++i) {
-        if ((in_u_dirichlet_arr[i])[node_lid]) {
-          DoFLocalId dof_lid = node_dof.dofId(node_lid, i);
-          Real u_dirichlet = in_u[node_lid][i];
-          /*if (ELIMINATE_BY_ROW_N_COLUMN)
-            linear_system.eliminateRowColumn(dof_lid, u_dirichlet);
-          else
-            linear_system.eliminateRow(dof_lid, u_dirichlet);*/
+    if constexpr (std::is_floating_point<DataType>::value) {
+      command << RUNCOMMAND_ENUMERATE(NodeLocalId, node_lid, m_mesh->ownNodes())
+      {
+        for (auto i = 0; i < NB_DOF; ++i) {
+          if ((in_u_dirichlet_arr[i])[node_lid]) {
+            DoFLocalId dof_lid = node_dof.dofId(node_lid, i);
+            Real u_dirichlet = in_u[node_lid];
+            if (!ELIMINATE_BY_ROW_N_COLUMN) {
+              in_out_dof_elimination_info[dof_lid] = 1; // How to define constant on device code ? ELIMINATE_ROW;
+              in_out_dof_elimination_value[dof_lid] = u_dirichlet;
+            }
+          }
         }
-      }
-    };
+      };
+    }
+    else
+      ARCANE_THROW(NotImplementedException, "BSRFormat(applyDirichletByElimination): Method not supported!");
+
+    linear_system->setEliminationArrays(m_dof_elimination_info, m_dof_elimination_value);
   }
 
   /*---------------------------------------------------------------------------*/
   /*---------------------------------------------------------------------------*/
 
-  void assembleLinearOperator(Arcane::String method, Real penalty, NumArray<Real, MDDim1>& rhs_vect, std::array<VariableNodeByte, NB_DOF>& u_dirichlet_arr, VariableNodeReal u)
+  template <typename DataType>
+  void applyDirichlet(Arcane::String method, Real penalty, NumArray<Real, MDDim1>& rhs_vect, std::array<VariableNodeByte, NB_DOF>& u_dirichlet_arr, MeshVariableScalarRefT<Node, DataType> u, DoFLinearSystem* linear_system)
   {
 
     if (method == "Penalty") {
-      m_use_csr_in_linear_system ? applyDirichletByPenalty<false, true>(penalty, rhs_vect, u_dirichlet_arr, u)
-                                 : applyDirichletByPenalty<false, false>(penalty, rhs_vect, u_dirichlet_arr, u);
+      m_use_csr_in_linear_system ? applyDirichletViaPenalty<false, true, DataType>(penalty, rhs_vect, u_dirichlet_arr, u)
+                                 : applyDirichletViaPenalty<false, false, DataType>(penalty, rhs_vect, u_dirichlet_arr, u);
     }
     else if (method == "WeaKPenalty") {
-      m_use_csr_in_linear_system ? applyDirichletByPenalty<true, true>(penalty, rhs_vect, u_dirichlet_arr, u)
-                                 : applyDirichletByPenalty<true, false>(penalty, rhs_vect, u_dirichlet_arr, u);
+      m_use_csr_in_linear_system ? applyDirichletViaPenalty<true, true, DataType>(penalty, rhs_vect, u_dirichlet_arr, u)
+                                 : applyDirichletViaPenalty<true, false, DataType>(penalty, rhs_vect, u_dirichlet_arr, u);
     }
+    else if (method == "RowElimination")
+      applyDirichletViaElimination<false, DataType>(u_dirichlet_arr, u, linear_system);
+    else if (method == "RowColumnElimination")
+      applyDirichletViaElimination<true, DataType>(u_dirichlet_arr, u, linear_system);
     else
       ARCANE_THROW(NotImplementedException, "BSRFormat(assembleLinearOperator): Method not supported!");
   }

femutils/DoFLinearSystem.cc

@@ -103,6 +103,11 @@ class SequentialDoFLinearSystemImpl
     ARCANE_THROW(NotImplementedException, "");
   }
 
+  void setEliminationArrays(VariableDoFByte& dof_elimination_info, VariableDoFReal& dof_elimination_value) override
+  {
+    ARCANE_THROW(NotImplementedException, "");
+  };
+
   void solve() override
   {
     _fillRHSVector();
@@ -379,6 +384,12 @@ eliminateRowColumn(DoFLocalId row, Real value)
   m_p->eliminateRowColumn(row, value);
 }
 
+void DoFLinearSystem::setEliminationArrays(VariableDoFByte& dof_elimination_info, VariableDoFReal& dof_elimination_value)
+{
+  _checkInit();
+  m_p->setEliminationArrays(dof_elimination_info, dof_elimination_value);
+};
+
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/
 

femutils/DoFLinearSystem.h

@@ -79,6 +79,7 @@ class DoFLinearSystemImpl
   virtual void matrixSetValue(DoFLocalId row, DoFLocalId column, Real value) = 0;
   virtual void eliminateRow(DoFLocalId row, Real value) = 0;
   virtual void eliminateRowColumn(DoFLocalId row, Real value) = 0;
+  virtual void setEliminationArrays(VariableDoFByte &dof_elimination_info, VariableDoFReal &dof_elimination_value) = 0;
   virtual void solve() = 0;
   virtual VariableDoFReal& solutionVariable() = 0;
   virtual VariableDoFReal& rhsVariable() = 0;
@@ -173,6 +174,8 @@ class DoFLinearSystem
    */
   void eliminateRowColumn(DoFLocalId row, Real value);
 
+  void setEliminationArrays(VariableDoFByte &dof_elimination_info, VariableDoFReal &dof_elimination_value);
+
   /*!
    * \brief Solve the current linear system.
    */

femutils/HypreDoFLinearSystem.cc

@@ -141,6 +141,11 @@ class HypreDoFLinearSystemImpl
     ARCANE_THROW(NotImplementedException, "");
   }
 
+  void setEliminationArrays(VariableDoFByte& dof_elimination_info, VariableDoFReal& dof_elimination_value) override
+  {
+    ARCANE_THROW(NotImplementedException, "");
+  };
+
   void solve() override;
 
   VariableDoFReal& solutionVariable() override

poisson/FemModule.cc

@@ -640,9 +640,20 @@ _assembleLinearOperator(BSRMatrix<1>* bsr_matrix)
   // Temporary variable to keep values for the RHS part of the linear system
   VariableDoFReal& rhs_values(m_linear_system.rhsVariable());
   rhs_values.fill(0.0);
-
   auto node_dof(m_dofs_on_nodes.nodeDoFConnectivityView());
 
+  if (m_use_csr) {
+  m_rhs_vect.resize(nbNode());
+  m_rhs_vect.fill(0.0);
+
+  std::array<VariableNodeByte, 1> u_dirichlet_arr = { m_u_dirichlet };
+
+  auto method = options()->enforceDirichletMethod();
+  m_bsr_format.applyDirichlet(options()->enforceDirichletMethod(), options()->penalty(), m_rhs_vect, u_dirichlet_arr, m_u, &m_linear_system);
+
+  _translateRhs();
+  }
+  else {
   if (options()->enforceDirichletMethod() == "Penalty") {
 
     Timer::Action timer_action(m_time_stats, "Penalty");
@@ -761,6 +772,7 @@ _assembleLinearOperator(BSRMatrix<1>* bsr_matrix)
            << "  - RowElimination\n"
            << "  - RowColumnElimination\n";
   }
+  }
 
   {
     Timer::Action timer_action(m_time_stats, "ConstantSourceTermAssembly");

poisson/FemModule.h

@@ -98,6 +98,7 @@
 #include <arcane/core/MeshUtils.h>
 #include "BSRFormat.h"
 #include "ArcaneFemFunctionsGpu.h"
+#include "arcane/core/VariableTypedef.h"
 
 /*---------------------------------------------------------------------------*/
 /*---------------------------------------------------------------------------*/