diff --git a/src/scheme/Order2AcousticSolver.cpp b/src/scheme/Order2AcousticSolver.cpp index c01a423cdca606474af757c357428cd02ebca264..719f80a5e81e699556b2ee04cfd7214493a86db0 100644 --- a/src/scheme/Order2AcousticSolver.cpp +++ b/src/scheme/Order2AcousticSolver.cpp @@ -631,7 +631,7 @@ class Order2AcousticSolverHandler::Order2AcousticSolver final : public Order2Aco const std::shared_ptr<const DiscreteFunctionVariant>& p_app = std::make_shared<const DiscreteFunctionVariant>(p_d); - const std::shared_ptr<const DiscreteFunctionVariant>& c_app = + const std::shared_ptr<const DiscreteFunctionVariant>& c_app = std::make_shared<const DiscreteFunctionVariant>(c_d); auto [ur_app, Fjr_app] = compute_fluxes(solver_type,rho_app,c_app,u_app,p_app,bc_descriptor_list); diff --git a/src/scheme/Order2Limiters.hpp b/src/scheme/Order2Limiters.hpp index 277e02d6b7116247c622850fd2bf50ae815b5f2d..74b2e2991ca84dac22e74a232c8a389d4dfc2669 100644 --- a/src/scheme/Order2Limiters.hpp +++ b/src/scheme/Order2Limiters.hpp @@ -87,6 +87,8 @@ scalar_limiter(const MeshType& mesh, coefficients[i] *= lambda; } }); + synchronize(DPk_fh.cellArrays()); + } template <MeshConcept MeshType> @@ -94,14 +96,14 @@ scalar_limiter(const MeshType& mesh, vector_limiter(const MeshType& mesh, const DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>& f, DiscreteFunctionDPk<MeshType::Dimension, TinyVector<MeshType::Dimension>> DPk_fh, - const CellValue<const TinyMatrix<MeshType::Dimension>>& grad_u) + const CellValue<const TinyMatrix<MeshType::Dimension>>& grad_u) { 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(); CellValue<TinyMatrix<MeshType::Dimension>> sym_grad_u{mesh.connectivity()}; @@ -170,7 +172,7 @@ scalar_limiter(const MeshType& mesh, } //} }); - + parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){ const double fn = dot(f[cell_id],n[cell_id]); @@ -227,12 +229,12 @@ scalar_limiter(const MeshType& mesh, 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); @@ -256,15 +258,16 @@ scalar_limiter(const MeshType& mesh, 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] + 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] + 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]; } }); + synchronize(DPk_fh.cellArrays()); } template <MeshConcept MeshType> @@ -285,7 +288,7 @@ scalar_limiter(const MeshType& mesh, parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id){ const double inv2j = inv2[cell_id]; - + double inv2_min = inv2j; double inv2_max = inv2j; @@ -311,7 +314,7 @@ scalar_limiter(const MeshType& mesh, double coef2 = 1.; if(std::abs(inv2_bar_max - inv2j) > eps){ coef1 = (inv2_max - inv2j) / (inv2_bar_max - inv2j); - } + } if(std::abs(inv2_bar_min - inv2j) > eps){ coef2 = (inv2_min - inv2j) / (inv2_bar_min - inv2j); } @@ -323,8 +326,9 @@ scalar_limiter(const MeshType& mesh, coefficients[0] = (1-lambda_inv2)*S[cell_id] + lambda_inv2 * coefficients[0]; for (size_t i = 1; i < coefficients.size(); ++i) { coefficients[i] *= lambda_inv2; - } + } }); - } + synchronize(DPk_Sh.cellArrays()); + } -#endif //ORDER2_LIMITERS_HPP \ No newline at end of file +#endif //ORDER2_LIMITERS_HPP