diff --git a/CMakeLists.txt b/CMakeLists.txt index 67136d4686b0caa3063ec9090af053ebd46ee136..e8e92ad5c72bbfeccb94c75a0dc10db9ffaa0d92 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -545,7 +545,6 @@ target_link_libraries( pugs PugsMesh PugsAlgebra - PugsAnalysis PugsUtils PugsLanguage PugsLanguageAST @@ -555,6 +554,7 @@ target_link_libraries( PugsAlgebra PugsAnalysis PugsScheme + PugsAnalysis PugsUtils PugsOutput PugsLanguageUtils diff --git a/src/analysis/Polynomial1D.hpp b/src/analysis/Polynomial1D.hpp index 894ffa8d8c03bb58f2389fc1b02917b181c759c5..ca19bcf9950d2fde179cf5f0b1fb82a8b3301185 100644 --- a/src/analysis/Polynomial1D.hpp +++ b/src/analysis/Polynomial1D.hpp @@ -263,6 +263,15 @@ class [[nodiscard]] Polynomial1D for (size_t i = 0; i < coeffients.size(); ++i) { m_coefficients[i] = coeffients[i]; } + + size_t real_degree = this->_getRealDegree(); + if (real_degree + 1 < m_coefficients.size()) { + Array<double> real_coefficients(real_degree + 1); + for (size_t i = 0; i <= real_degree; ++i) { + real_coefficients[i] = m_coefficients[i]; + } + m_coefficients = real_coefficients; + } } PUGS_INLINE diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index b4022d660475cf93eb9f9e1cb291622de22dba22..c6036fda9e1873b285a4a40fd96b2833bfadc64c 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -188,6 +188,42 @@ SchemeModule::SchemeModule() )); + this->_addBuiltinFunction("entropic_eucclhyd_solver", + std::make_shared<BuiltinFunctionEmbedder<std::tuple< + double, std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&, + const std::vector<std::shared_ptr< + const IBoundaryConditionDescriptor>>&, + const double&)>>( + + [](const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& E, + const std::shared_ptr<const IDiscreteFunction>& c, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& + bc_descriptor_list, + const double& dt_max) + -> std::tuple< + double, std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, std::shared_ptr<const IDiscreteFunction>> { + return AcousticSolverHandler{AcousticSolverHandler::SolverType::Eucclhyd, + rho, + c, + u, + p, + bc_descriptor_list} + .solver() + .apply_computing_entropic_dt(dt_max, rho, u, E, p); + } + + )); + this->_addBuiltinFunction("eucclhyd_solver", std::make_shared<BuiltinFunctionEmbedder<std::tuple< std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>, diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp index 6973037571b23f861cec5c70ed28fbe67c036fcf..1d04fb5a5e8a12ab5444e7421ffb61f2649e8261 100644 --- a/src/scheme/AcousticSolver.cpp +++ b/src/scheme/AcousticSolver.cpp @@ -1,5 +1,7 @@ #include <scheme/AcousticSolver.hpp> +#include <analysis/Polynomial1D.hpp> +#include <analysis/Polynomial1DRootComputer.hpp> #include <language/utils/InterpolateItemValue.hpp> #include <mesh/ItemValueUtils.hpp> #include <mesh/MeshNodeBoundary.hpp> @@ -84,6 +86,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler NodeValue<const Rd> m_ur; NodeValuePerCell<const Rd> m_Fjr; + NodeValuePerCell<const Rdxd> m_Ajr; + CellValue<const double> _getRhoC(const DiscreteScalarFunction& rho, const DiscreteScalarFunction& c) { @@ -393,6 +397,294 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler return F; } + double + _computeMaxDensityDt(const MeshType& mesh, + const DiscreteFunctionP0<Dimension, double>& rho, + const double dt_max) const + { + double dt = dt_max; + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + const CellValue<const double> Vj = mesh_data.Vj(); + + const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + + CellValue<double> dVj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double dV = 0; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + dV += dot(Cjr(j, R), m_ur[r]); + } + + dVj[j] = dV; + }); + + CellValue<double> Gj{mesh.connectivity()}; + + if constexpr (Dimension == 2) { + for (CellId j = 0; j < mesh.numberOfCells(); ++j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double G = 0; + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()]; + G -= dot(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]}); + } + Gj[j] = 0.5 * G; + } + } + + for (CellId j = 0; j < mesh.numberOfCells(); ++j) { + Polynomial1D p_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])}); + if constexpr (Dimension == 2) { + p_tau = p_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])}); + } + Polynomial1DRootComputer computer{p_tau, 0, dt}; + std::optional dt_tau = computer.getFirstRoot(); + if (dt_tau.has_value()) { + Assert(dt_tau.value() <= dt); + dt = dt_tau.value(); + } + } + + if (dt != dt_max) { + std::cout << "volume variation imposes time step\n"; + } + + if (dt < dt_max) { + dt *= 0.95 * dt; + } + + return parallel::allReduceMin(dt); + } + + double + _computeMaxEnergyDt(const MeshType& mesh, + const DiscreteScalarFunction& rho, + const DiscreteVectorFunction& u, + const DiscreteScalarFunction& E, + const DiscreteScalarFunction& p, + const double dt_max) const + { + double dt = dt_max; + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + const CellValue<const double> Vj = mesh_data.Vj(); + + const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + + CellValue<double> dVj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double dV = 0; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + dV += dot(Cjr(j, R), m_ur[r]); + } + + dVj[j] = dV; + }); + + CellValue<double> dSj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double dS = 0; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + + Rd du = u[j] - m_ur[r]; + dS += dot(m_Ajr(j, R) * du, du); + } + + dSj[j] = dS; + }); + + CellValue<Rd> dPj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + Rd dP = zero; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + + Rd du = u[j] - m_ur[r]; + dP += m_Ajr(j, R) * du; + } + + dPj[j] = dP; + }); + + for (CellId j = 0; j < mesh.numberOfCells(); ++j) { + const double inv_mj = 1 / (rho[j] * Vj[j]); + + Polynomial1D p_e( + {E[j] - 0.5 * dot(u[j], u[j]), inv_mj * (dSj[j] - p[j] * dVj[j]), -inv_mj * inv_mj * dot(dPj[j], dPj[j])}); + Polynomial1DRootComputer computer{p_e, 0, dt}; + std::optional dt_e = computer.getFirstRoot(); + if (dt_e.has_value()) { + Assert(dt_e.value() <= dt); + dt = dt_e.value(); + } + } + + if (dt < dt_max) { + dt *= 0.95 * dt; + } + + if (dt != dt_max) { + std::cout << "internal energy variation imposes time step\n"; + } + + return parallel::allReduceMin(dt); + } + + double + _computeMaxEntropyDt(const MeshType& mesh, + const DiscreteScalarFunction& rho, + const DiscreteVectorFunction& u, + const DiscreteScalarFunction& p, + const double dt_max) const + { + double dt = dt_max; + + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + const CellValue<const double> Vj = mesh_data.Vj(); + + const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + + CellValue<double> dVj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double dV = 0; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + dV += dot(Cjr(j, R), m_ur[r]); + } + + dVj[j] = dV; + }); + + CellValue<double> Gj{mesh.connectivity()}; + + if constexpr (Dimension == 2) { + for (CellId j = 0; j < mesh.numberOfCells(); ++j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double G = 0; + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()]; + G -= dot(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]}); + } + Gj[j] = 0.5 * G; + } + } + + CellValue<double> dSj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + double dS = 0; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + + const Rd du = u[j] - m_ur[r]; + dS += dot(m_Ajr(j, R) * du, du); + } + + dSj[j] = dS; + }); + + CellValue<Rd> dPj{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + Rd dP = zero; + + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + + Rd du = u[j] - m_ur[r]; + dP += m_Ajr(j, R) * du; + } + + dPj[j] = dP; + }); + +#warning fixed gamma value + const double gamma = 1.4; + + for (CellId j = 0; j < mesh.numberOfCells(); ++j) { + const double inv_mj = 1 / (rho[j] * Vj[j]); + const double inv_Vj = 1 / Vj[j]; + + Polynomial1D DT({0, 1}); + std::vector<double> delta_S = {dSj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])}; + if constexpr (Dimension == 2) { + delta_S[1] += p[j] * Gj[j]; + } + Polynomial1D p_S0(delta_S); + + std::vector<double> delta_V{dVj[j] /*, dG */}; + if constexpr (Dimension == 2) { + delta_V.push_back(Gj[j]); + } + Polynomial1D p_deltaV(delta_V); + + Polynomial1D p_S1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])}); + + Polynomial1D p_S = // + p_S0 // + + (inv_Vj * DT) * (((gamma - 1) * p_S1) * p_deltaV + (gamma - 2) * p[j] * p_deltaV * p_deltaV) // + + (inv_Vj * DT) * (inv_Vj * DT) * (((gamma - 1) * (gamma - 2)) * p_S1) * p_deltaV * p_deltaV; + + Polynomial1DRootComputer computer{p_S, 0, dt}; + std::optional dt_S = computer.getFirstRoot(); + if (dt_S.has_value()) { + Assert(dt_S.value() <= dt); + dt = dt_S.value(); + } + } + + if (dt != dt_max) { + std::cout << "entropy variation imposes time step\n"; + } + return parallel::allReduceMin(dt); + } + AcousticSolver(SolverType solver_type, const std::shared_ptr<const MeshType>& p_mesh, const DiscreteScalarFunction& rho, @@ -407,6 +699,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler CellValue<const double> rhoc = this->_getRhoC(rho, c); NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rhoc); + m_Ajr = Ajr; + NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr); NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p); @@ -494,6 +788,91 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E)); } + std::tuple<double, + std::shared_ptr<const IMesh>, + std::shared_ptr<const DiscreteScalarFunction>, + std::shared_ptr<const DiscreteVectorFunction>, + std::shared_ptr<const DiscreteScalarFunction>> + apply_computing_entropic_dt(const double& dt_max, + const std::shared_ptr<const MeshType>& mesh, + const DiscreteScalarFunction& rho, + const DiscreteVectorFunction& u, + const DiscreteScalarFunction& E, + const DiscreteScalarFunction& p) const + { + const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + + double dt = dt_max; + + dt = _computeMaxDensityDt(*mesh, rho, dt); + dt = _computeMaxEnergyDt(*mesh, rho, u, E, p, dt); + dt = _computeMaxEntropyDt(*mesh, rho, u, p, dt); + + NodeValue<Rd> new_xr = copy(mesh->xr()); + parallel_for( + mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; }); + + std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr); + + CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj(); + + CellValue<double> new_rho = copy(rho.cellValues()); + CellValue<Rd> new_u = copy(u.cellValues()); + CellValue<double> new_E = copy(E.cellValues()); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + Rd momentum_fluxes = zero; + double energy_fluxes = 0; + for (size_t R = 0; R < cell_nodes.size(); ++R) { + const NodeId r = cell_nodes[R]; + momentum_fluxes += m_Fjr(j, R); + energy_fluxes += dot(m_Fjr(j, R), m_ur[r]); + } + const double dt_over_Mj = dt / (rho[j] * Vj[j]); + new_u[j] -= dt_over_Mj * momentum_fluxes; + new_E[j] -= dt_over_Mj * energy_fluxes; + }); + + CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; }); + + return {dt, new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho), + std::make_shared<DiscreteVectorFunction>(new_mesh, new_u), + std::make_shared<DiscreteScalarFunction>(new_mesh, new_E)}; + } + + std::tuple<double, + std::shared_ptr<const IMesh>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>> + apply_computing_entropic_dt(const double& dt_max, + const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& E, + const std::shared_ptr<const IDiscreteFunction>& p) const + { + std::shared_ptr i_mesh = getCommonMesh({rho, u, E}); + if (not i_mesh) { + throw NormalError("discrete functions are not defined on the same mesh"); + } + + if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) { + throw NormalError("acoustic solver expects P0 functions"); + } + + return this->apply_computing_entropic_dt(dt_max, std::dynamic_pointer_cast<const MeshType>(i_mesh), + *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho), + *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u), + *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E), + *std::dynamic_pointer_cast<const DiscreteScalarFunction>(p)); + } + AcousticSolver(SolverType solver_type, const std::shared_ptr<const IMesh>& mesh, const std::shared_ptr<const IDiscreteFunction>& rho, diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 2ed97075717f8984afb22799648b4b26fbc9327b..8866e09f8c740895e5dccf88a5b5ceed0de2ce88 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -32,6 +32,17 @@ class AcousticSolverHandler const std::shared_ptr<const IDiscreteFunction>& u, const std::shared_ptr<const IDiscreteFunction>& E) const = 0; + virtual std::tuple<double, + std::shared_ptr<const IMesh>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>> + apply_computing_entropic_dt(const double& dt_max, + const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& E, + const std::shared_ptr<const IDiscreteFunction>& p) const = 0; + IAcousticSolver() = default; IAcousticSolver(IAcousticSolver&&) = default; IAcousticSolver& operator=(IAcousticSolver&&) = default; diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt index 2131f2fec0748407ea7bb21a3375d398e1f6242a..3a346f30dc129e7ba76e1d503b8d227a05906b5f 100644 --- a/src/scheme/CMakeLists.txt +++ b/src/scheme/CMakeLists.txt @@ -6,3 +6,9 @@ add_library( DiscreteFunctionInterpoler.cpp DiscreteFunctionUtils.cpp DiscreteFunctionVectorInterpoler.cpp) + +# Additional dependencies +add_dependencies( + PugsScheme + PugsAnalysis +) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4a602a294bf663d41e075faf29ebaa943058ca54..2e75a839324baf50f9903cfb16c528f4cab384be 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -130,8 +130,8 @@ target_link_libraries (unit_tests PugsLanguage PugsMesh PugsAlgebra - PugsAnalysis PugsScheme + PugsAnalysis PugsOutput PugsUtils kokkos @@ -153,10 +153,10 @@ target_link_libraries (mpi_unit_tests PugsLanguageAlgorithms PugsMesh PugsAlgebra - PugsAnalysis PugsUtils PugsLanguageUtils PugsScheme + PugsAnalysis PugsOutput PugsUtils PugsAlgebra