Add new velocity limiter

src/scheme/Order2AcousticSolver.cpp

@@ -409,73 +409,175 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
   void 
   _vector_limiter(const MeshType& mesh,
                   const DiscreteFunctionP0<const Rd>& f,
-                  DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) const
+                  DiscreteFunctionDPk<Dimension, Rd>& DPK_fh,
+                  //const DiscreteFunctionP0<const double>& p,
+                  const DiscreteFunctionDPk<Dimension, const double>& DPk_ph) const 
   {
     MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
     StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
     StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
-    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor, symmetry_boundary_descriptor_list);
+    auto stencil = StencilManager::instance().getCellToCellStencilArray(mesh.connectivity(), stencil_descriptor,
+                                                                        symmetry_boundary_descriptor_list);
     
     const auto xj = mesh_data.xj();
 
-    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
-      const Rd fj = f[cell_id];
+    //Calcul du gradiant
+    CellValue<Rd> n{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        auto coefficients_p = DPk_ph.coefficients(cell_id);
+        Rd grad_p = zero;
+        for(size_t i = 1; i < coefficients_p.size(); ++i){
+          grad_p[i] = coefficients_p[i];
+        }
+        const double norm_grad_p = l2Norm(grad_p);
+        if(norm_grad_p == 0){
+          n[cell_id] = zero;
+        }
+        else{
+          for(size_t d = 0; d < Dimension; ++d){
+            n[cell_id][d] = grad_p[d] / norm_grad_p;
+          }
+        }
 
-      Rd f_min = fj;
-      Rd f_max = fj;
+        
+      });
+
+
+    //Construction des vecteurs orthogonaux
+    const CellValue<Rd> t{mesh.connectivity()};
+    const CellValue<Rd> l{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const Rd nj = n[cell_id];
+        Rd a = zero;
+        if(l2Norm(nj) != 0){
+          if((nj[0] / l2Norm(nj) != 1) and (nj[0] / l2Norm(nj) != -1)){
+            a[0] = 1.;
+          }
+          else { 
+            a[1] = 1.;
+          }
+        }
+
+        Rd tj = a - dot(a,nj) * nj;
+        const double& norm_tj = l2Norm(tj);
+        if(norm_tj == 0){
+          tj = zero;
+        }
+        else {
+          for(size_t d = 0; d < Dimension; ++d){
+            tj[d] = tj[d] / norm_tj;
+          }
+        }
+        t[cell_id] = tj;
+
+        Rd lj = zero;
+        if(Dimension == 3){
+          lj[0] = nj[1]*tj[2] - nj[2]*tj[1];
+          lj[1] = nj[2]*tj[0] - nj[0]*tj[2];
+          lj[2] = nj[0]*tj[1] - nj[1]*tj[0];
+          const double& norm_lj = l2Norm(lj);
+          for(size_t d = 0; d < Dimension; ++d){
+            lj[d] = lj[d] / norm_lj;
+          }
+        }
+        l[cell_id] = lj;
+      });
+
+    //Limiteurs
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+        const double fn = dot(f[cell_id],n[cell_id]);
+        const double ft = dot(f[cell_id],t[cell_id]);
+        const double fl = dot(f[cell_id],l[cell_id]);
+
+        double fn_min = fn;
+        double fn_max = fn;
+        double ft_min = ft;
+        double ft_max = ft;
+        double fl_min = fl;
+        double fl_max = fl;
 
         const auto cell_stencil = stencil[cell_id];
         for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
-        for(size_t dim = 0; dim < Dimension; ++dim){
-          f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][dim]);
-          f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][dim]);
-        }
+          const double fn_k = dot(f[cell_stencil[i_cell]],n[cell_stencil[i_cell]]);
+          fn_min            = std::min(fn_min,fn_k);
+          fn_max            = std::max(fn_max,fn_k);
+
+          const double ft_k = dot(f[cell_stencil[i_cell]],t[cell_stencil[i_cell]]);
+          ft_min            = std::min(ft_min,ft_k);
+          ft_max            = std::max(ft_max,ft_k);
+
+          const double fl_k = dot(f[cell_stencil[i_cell]],l[cell_stencil[i_cell]]);
+          fl_min            = std::min(fl_min,fl_k);
+          fl_max            = std::max(fl_max,fl_k);
         }
 
-      Rd f_bar_min = fj;
-      Rd f_bar_max = fj;
+        double fn_bar_min = fn;
+        double fn_bar_max = fn;
+        double ft_bar_min = ft;
+        double ft_bar_max = ft;
+        double fl_bar_min = fl;
+        double fl_bar_max = fl;
 
         for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
           const CellId cell_k_id = cell_stencil[i_cell];
-        const Rd f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
+          const Rd f_xk          = DPK_fh[cell_id](xj[cell_k_id]);
 
-        for(size_t dim = 0; dim < Dimension; ++dim){
-          f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]);
-          f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]);
-        }
+          const double fn_xk = dot(f_xk,n[cell_k_id]);
+          fn_bar_min         = std::min(fn_bar_min,fn_xk);
+          fn_bar_max         = std::max(fn_bar_max,fn_xk);
+
+          const double ft_xk = dot(f_xk,t[cell_k_id]);
+          ft_bar_min         = std::min(ft_bar_min,ft_xk);
+          ft_bar_max         = std::max(ft_bar_max,ft_xk);
+
+          const double fl_xk = dot(f_xk,l[cell_k_id]);
+          fl_bar_min         = std::min(fl_bar_min,fl_xk);
+          fl_bar_max         = std::max(fl_bar_max,fl_xk);
         }
 
         const double eps = 1E-14;
-      Rd coef1;
-      for(size_t dim = 0; dim < Dimension; ++dim){
-        coef1[dim] = 1;
-        if (std::abs(f_bar_max[dim] - fj[dim]) > eps) {
-          coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim]));
+        double coef1_n   = 1;
+        if(std::abs(fn_bar_max - fn) > eps){
+          coef1_n = (fn_max - fn) / (fn_bar_max - fn);
         } 
+        double coef2_n = 1;
+        if(std::abs(fn_bar_min - fn) > eps){
+          coef2_n = (fn_min - fn) / (fn_bar_min - fn);
         }
         
-      Rd coef2;
-      for(size_t dim = 0; dim < Dimension; ++dim){
-        coef2[dim] = 1;
-        if (std::abs(f_bar_min[dim] - fj[dim]) > eps) {
-          coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim]));
+        double coef1_t = 1;
+        if(std::abs(ft_bar_max - ft) > eps){
+          coef1_t = (ft_max - ft) / (ft_bar_max - ft);
         }
+        double coef2_t = 1;
+        if(std::abs(ft_bar_min - ft) > eps){
+          coef2_t = (ft_min - ft) / (ft_bar_min - ft);
         }
 
-      double min_coef1 = coef1[0];
-      double min_coef2 = coef2[0];
-      for(size_t dim = 1; dim < Dimension; ++dim){
-        min_coef1 = std::min(min_coef1,coef1[dim]);
-        min_coef2 = std::min(min_coef2,coef2[dim]);
+        double coef1_l = 1;
+        if(std::abs(fl_bar_max - fl) > eps){
+          coef1_l = (fl_max - fl) / (fl_bar_max - fl);
+        }
+        double coef2_l = 1;
+        if(std::abs(fl_bar_min - fl) > eps){
+          coef2_l = (fl_min - fl) / (fl_bar_min - fl);
         }
 
-      const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
-
-      auto coefficients = DPk_fh.coefficients(cell_id);
+        const double lambda_n = std::max(0., std::min(1., std::min(coef1_n,coef2_n)));
+        const double lambda_t = std::max(0., std::min(1., std::min(coef1_t,coef2_t)));
+        const double lambda_l = std::max(0., std::min(1., std::min(coef1_l,coef2_l)));
 
-      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+        auto coefficients = DPK_fh.coefficients(cell_id);
+        coefficients[0] = (1 - lambda_n)*fn*n[cell_id] + lambda_n*dot(coefficients[0],n[cell_id])*n[cell_id] 
+                        + (1 - lambda_t)*ft*t[cell_id] + lambda_t*dot(coefficients[0],t[cell_id])*t[cell_id] 
+                        + (1 - lambda_l)*fl*l[cell_id] + lambda_l*dot(coefficients[0],l[cell_id])*l[cell_id];
         for(size_t i = 1; i < coefficients.size(); ++i){
-        coefficients[i] *= lambda;
+          coefficients[i] = lambda_n*dot(coefficients[i],n[cell_id])*n[cell_id] 
+                          + lambda_t*dot(coefficients[i],t[cell_id])*t[cell_id] 
+                          + lambda_l*dot(coefficients[i],l[cell_id])*l[cell_id];
         }
     });
   }
@@ -527,7 +629,7 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_ph);
 
     this->_scalar_limiter(mesh,p,p_lim);
-    this->_vector_limiter(mesh,u,u_lim);
+    this->_vector_limiter(mesh,u,u_lim,DPk_ph);
 
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);