Improve periodic boundary conditions treatment

src/scheme/EulerKineticSolverLagrangeMultiDLocalOrder2Periodic.cpp

@@ -61,6 +61,8 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
     const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
     const auto& node_to_cell_matrix               = p_mesh->connectivity().nodeToCellMatrix();
 
+    NodeArray<double> sum_j_li_njr(m_mesh.connectivity(), nb_waves);
+
     Fr.fill(0.);
 
     parallel_for(
@@ -97,6 +99,26 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
               sum_li_njr += li_njr;
             }
           }
+          sum_j_li_njr[node_id][i_wave] = sum_li_njr;
+        }
+      });
+
+    // Periodic BC
+    for (size_t i = 0; i < Fr.sizeOfArrays(); ++i) {
+      auto F                                    = Fr[NodeId{0}][i] + Fr[NodeId(m_mesh.numberOfNodes() - 1)][i];
+      Fr[NodeId{0}][i]                          = F;
+      Fr[NodeId(m_mesh.numberOfNodes() - 1)][i] = F;
+
+      auto s = sum_j_li_njr[NodeId{0}][i] + sum_j_li_njr[NodeId(m_mesh.numberOfNodes() - 1)][i];
+
+      sum_j_li_njr[NodeId{0}][i]                          = s;
+      sum_j_li_njr[NodeId(m_mesh.numberOfNodes() - 1)][i] = s;
+    }
+
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+          const double sum_li_njr = sum_j_li_njr[node_id][i_wave];
           if (sum_li_njr != 0) {
             Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
           }
@@ -117,6 +139,8 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
     const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
     const auto& node_to_cell_matrix               = p_mesh->connectivity().nodeToCellMatrix();
 
+    NodeArray<double> sum_j_li_njr(m_mesh.connectivity(), nb_waves);
+
     Fr.fill(zero);
 
     parallel_for(
@@ -151,6 +175,26 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
               sum_li_njr += li_njr;
             }
           }
+          sum_j_li_njr[node_id][i_wave] = sum_li_njr;
+        }
+      });
+
+    // Periodic BC
+    for (size_t i = 0; i < Fr.sizeOfArrays(); ++i) {
+      auto F                                    = Fr[NodeId{0}][i] + Fr[NodeId(m_mesh.numberOfNodes() - 1)][i];
+      Fr[NodeId{0}][i]                          = F;
+      Fr[NodeId(m_mesh.numberOfNodes() - 1)][i] = F;
+
+      auto s = sum_j_li_njr[NodeId{0}][i] + sum_j_li_njr[NodeId(m_mesh.numberOfNodes() - 1)][i];
+
+      sum_j_li_njr[NodeId{0}][i]                          = s;
+      sum_j_li_njr[NodeId(m_mesh.numberOfNodes() - 1)][i] = s;
+    }
+
+    parallel_for(
+      m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+          const double sum_li_njr = sum_j_li_njr[node_id][i_wave];
           if (sum_li_njr != 0) {
             Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
           }
@@ -176,6 +220,7 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
         }
         u[node_id] = inv_lambda_squared * u[node_id];
       });
+
     return u;
   }
 
@@ -322,6 +367,13 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
         }
 
         if constexpr (Dimension == 1) {
+          if (node_id == 0) {
+            lambda_r = std::max(lambda_r, lambda_j[CellId(m_mesh.numberOfCells() - 1)]);
+          }
+          if (node_id == m_mesh.numberOfNodes() - 1) {
+            lambda_r = std::max(lambda_r, lambda_j[CellId(0)]);
+          }
+
           for (size_t i_wave = 0; i_wave < m_k; ++i_wave) {
             lambda_vector[node_id][i_wave][0] = lambda_r * (1 - 2. * i_wave / (m_k - 1));
           }
@@ -349,7 +401,6 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
       });
     return lambda_vector;
   }
-
   void
   scalar_limiter(const CellValue<const double>& u, DiscreteFunctionDPk<Dimension, double>& DPk_uh) const
   {
@@ -451,7 +502,7 @@ class EulerKineticSolverLagrangeMultiDLocalOrder2Periodic
              std::shared_ptr<const DiscreteFunctionVariant>>
   apply()
   {
-    const bool limitation = true;
+    const bool limitation = false;
 
     const CellValue<double> rho              = copy(m_rho);
     const CellValue<TinyVector<Dimension>> u = copy(m_u);