From 7fd7d9614b3ac267367f0ab9c8dcb67d0c1ef949 Mon Sep 17 00:00:00 2001
From: Compte local pour Alexandre GANGLOFF <alexandre.gangloff@cea.fr>
Date: Wed, 6 Nov 2024 17:57:12 +0100
Subject: [PATCH] Add scalar and vector limiters for AcousticSolveur
---
src/scheme/Order2AcousticSolver.cpp | 143 +++++++++++++++++++++++++++-
1 file changed, 140 insertions(+), 3 deletions(-)
diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp
index d59676359..92f13f2c0 100644
--- a/src/scheme/Order2AcousticSolver.cpp
+++ b/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>();
--
GitLab