diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp
index 92f13f2c07c0c2af9ce0b7712fa4a41a5f40068d..c8f2d9e7bde77263421e5999584c16195a9e1afb 100644
--- a/src/scheme/Order2AcousticSolver.cpp
+++ b/src/scheme/Order2AcousticSolver.cpp
@@ -406,77 +406,179 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco
});
}
- void
+ 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];
-
- 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]);
+ //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_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]);
+ //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.;
+ }
}
- }
- 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 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;
+ });
- 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]));
+ //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){
+ 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);
}
- }
- 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 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;
- const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2)));
+ 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]);
- auto coefficients = DPk_fh.coefficients(cell_id);
+ 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);
- coefficients[0] = (1 - lambda) * f[cell_id] + lambda * coefficients[0];
- for (size_t i = 1; i < coefficients.size(); ++i) {
- coefficients[i] *= lambda;
- }
+ 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;
+ 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);
+ }
+
+ 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 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_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)));
+
+ 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_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);