Add scalar and vector limiters for AcousticSolveur

src/scheme/Order2AcousticSolver.cpp

@@ -349,6 +349,137 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     return F;
   }
 
+  void
+  _scalar_limiter(const MeshType& mesh,
+                  const DiscreteFunctionP0<const double>& f,
+                  DiscreteFunctionDPk<Dimension, double>& DPk_fh) 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);
+
+    const auto xj = mesh_data.xj();
+
+    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+      const double fj = f[cell_id];
+
+      double f_min = fj;
+      double f_max = fj;
+
+      const auto cell_stencil = stencil[cell_id];
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        f_min = std::min(f_min, f[cell_stencil[i_cell]]);
+        f_max = std::max(f_max, f[cell_stencil[i_cell]]);
+      }
+
+      double f_bar_min = fj;
+      double f_bar_max = fj;
+
+      for(size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell){
+        const CellId cell_k_id = cell_stencil[i_cell];
+        const double f_xk    = DPk_fh[cell_id](xj[cell_k_id]);
+
+        f_bar_min = std::min(f_bar_min, f_xk);
+        f_bar_max = std::max(f_bar_max, f_xk);
+      }
+
+      const double eps = 1E-14;
+      double coef1     = 1;
+      if (std::abs(f_bar_max - fj) > eps) {
+        coef1 = (f_max - fj) / ((f_bar_max - fj));
+      }
+
+      double coef2 = 1.;
+      if (std::abs(f_bar_min - fj) > eps) {
+        coef2 = (f_min - fj) / ((f_bar_min - fj));
+      }
+
+      const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+      auto coefficients = DPk_fh.coefficients(cell_id);
+
+      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+      for (size_t i = 1; i < coefficients.size(); ++i) {
+        coefficients[i] *= lambda;
+      }
+    });
+  }
+
+  void
+  _vector_limiter(const MeshType& mesh,
+                  const DiscreteFunctionP0<const Rd>& f,
+                  DiscreteFunctionDPk<Dimension, Rd>& DPk_fh) 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);
+
+    const auto xj = mesh_data.xj();
+
+    parallel_for(mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){
+      const Rd fj = f[cell_id];
+
+      Rd f_min = fj;
+      Rd f_max = fj;
+
+      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]);
+        }
+      }
+
+      Rd f_bar_min = fj;
+      Rd f_bar_max = fj;
+
+      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]);
+
+        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 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]));
+        }
+      }
+
+      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 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]);
+      }
+
+      const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
+
+      auto coefficients = DPk_fh.coefficients(cell_id);
+
+      coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
+      for (size_t i = 1; i < coefficients.size(); ++i) {
+        coefficients[i] *= lambda;
+      }
+    });
+  }
+
  public:
   std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
@@ -392,10 +523,16 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     auto DPk_uh = reconstruction[0]->get<DiscreteFunctionDPk<Dimension, const Rd>>();
     auto DPk_ph = reconstruction[1]->get<DiscreteFunctionDPk<Dimension, const double>>();
 
+    DiscreteFunctionDPk<Dimension, Rd> u_lim     = copy(DPk_uh);
+    DiscreteFunctionDPk<Dimension, double> p_lim = copy(DPk_ph);
+
+    this->_scalar_limiter(mesh,p,p_lim);
+    this->_vector_limiter(mesh,u,u_lim);
+
     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, DPk_uh, DPk_ph);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u_lim, p_lim);
 
     const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
     this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
@@ -405,7 +542,7 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     synchronize(br);
 
     NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
-    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, DPk_uh, DPk_ph);
+    NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u_lim, p_lim);
 
     return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -519,7 +656,7 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
     auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
     std::tie(m_app,rho_app,u_app,E_app) = apply_fluxes(dt/2, rho, u, E, ur, Fjr);
 
-    const double& gamma = 3;
+    const double& gamma = 1.4;
 
     const DiscreteScalarFunction& rho_d = rho_app->get<DiscreteScalarFunction>();
     const DiscreteVectorFunction& u_d   = u_app->get<DiscreteVectorFunction>();