diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp index edce9e025f0feb4cd39acaadcb537451c9251065..0af71d5a5eb7e07db2625bfa405309a47213fcf2 100644 --- a/src/language/algorithms/AcousticSolverAlgorithm.cpp +++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp @@ -24,7 +24,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm( std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); - typename AcousticSolver<MeshType>::BoundaryConditionList bc_list; + typename AcousticSolverOld<MeshType>::BoundaryConditionList bc_list; { constexpr ItemType FaceType = [] { if constexpr (Dimension > 1) { @@ -39,7 +39,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm( switch (bc_descriptor->type()) { case IBoundaryConditionDescriptor::Type::symmetry: { - using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition; + using SymmetryBoundaryCondition = typename AcousticSolverOld<MeshType>::SymmetryBoundaryCondition; const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor); @@ -60,7 +60,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm( const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor); if (dirichlet_bc_descriptor.name() == "velocity") { - using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition; + using VelocityBoundaryCondition = typename AcousticSolverOld<MeshType>::VelocityBoundaryCondition; for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { @@ -80,7 +80,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm( } } } else if (dirichlet_bc_descriptor.name() == "pressure") { - using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition; + using PressureBoundaryCondition = typename AcousticSolverOld<MeshType>::PressureBoundaryCondition; for (size_t i_ref_face_list = 0; i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { @@ -141,7 +141,7 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm( } unknowns.gammaj().fill(1.4); - AcousticSolver acoustic_solver(mesh, bc_list, solver_type); + AcousticSolverOld acoustic_solver(mesh, bc_list, solver_type); const double tmax = 0.2; double t = 0; diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index c23d3794e4c26a95301599b6933fe0d643391010..212464f79ce7ccdf24c0164aab84fe5c9b527c31 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -4,9 +4,11 @@ #include <language/utils/BuiltinFunctionEmbedder.hpp> #include <language/utils/TypeDescriptor.hpp> #include <mesh/Mesh.hpp> +#include <scheme/AcousticSolver.hpp> #include <scheme/DirichletBoundaryConditionDescriptor.hpp> #include <scheme/DiscreteFunctionDescriptorP0.hpp> #include <scheme/DiscreteFunctionInterpoler.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> #include <scheme/FourierBoundaryConditionDescriptor.hpp> #include <scheme/FreeBoundaryConditionDescriptor.hpp> #include <scheme/IBoundaryConditionDescriptor.hpp> @@ -18,6 +20,8 @@ #include <scheme/NumberedBoundaryDescriptor.hpp> #include <scheme/SymmetryBoundaryConditionDescriptor.hpp> +#include <scheme/PerfectGas.hpp> + #include <memory> SchemeModule::SchemeModule() @@ -201,5 +205,127 @@ SchemeModule::SchemeModule() } } + )); + + this + ->_addBuiltinFunction("perfect_gas_epsilon_from_rho_p_gamma", + std::make_shared<BuiltinFunctionEmbedder< + 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>& rho, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::shared_ptr<const IDiscreteFunction>& gamma) + -> std::shared_ptr<const IDiscreteFunction> { + return perfect_gas::epsilonFromRhoPGamma(rho, p, gamma); + } + + )); + + this + ->_addBuiltinFunction("perfect_gas_p_from_rho_epsilon_gamma", + std::make_shared<BuiltinFunctionEmbedder< + 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>& rho, + const std::shared_ptr<const IDiscreteFunction>& epsilon, + const std::shared_ptr<const IDiscreteFunction>& gamma) + -> std::shared_ptr<const IDiscreteFunction> { + return perfect_gas::pFromRhoEpsilonGamma(rho, epsilon, gamma); + } + + )); + + this + ->_addBuiltinFunction("perfect_gas_sound_speed", + std::make_shared<BuiltinFunctionEmbedder< + 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>& rho, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::shared_ptr<const IDiscreteFunction>& gamma) + -> std::shared_ptr<const IDiscreteFunction> { + return perfect_gas::soundSpeed(rho, p, gamma); + } + + )); + + this + ->_addBuiltinFunction("E_from_epsilon_u", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&)>>( + + [](const std::shared_ptr<const IDiscreteFunction>& epsilon, + const std::shared_ptr<const IDiscreteFunction>& u) + -> std::shared_ptr<const IDiscreteFunction> { + return perfect_gas::totalEnergyFromEpsilonU(epsilon, u); + } + + )); + + this + ->_addBuiltinFunction("epsilon_from_E_u", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&, + const std::shared_ptr<const IDiscreteFunction>&)>>( + + [](const std::shared_ptr<const IDiscreteFunction>& E, + const std::shared_ptr<const IDiscreteFunction>& u) + -> std::shared_ptr<const IDiscreteFunction> { + return perfect_gas::epsilonFromTotalEnergyU(E, u); + } + + )); + + this->_addBuiltinFunction( + "glace_solver", + std::make_shared<BuiltinFunctionEmbedder<std::tuple< + 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) + -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>> { + return AcousticSolverHandler{rho, c, u, p, bc_descriptor_list}.solver().apply(dt, rho, u, E); + } + + )); + + this + ->_addBuiltinFunction("lagrangian", + std::make_shared<BuiltinFunctionEmbedder< + std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IMesh>&, + const std::shared_ptr<const IDiscreteFunction>&)>>( + + [](const std::shared_ptr<const IMesh>& mesh, + const std::shared_ptr<const IDiscreteFunction>& v) + -> std::shared_ptr<const IDiscreteFunction> { return shallowCopy(mesh, v); } + + )); + + this->_addBuiltinFunction("acoustic_dt", + std::make_shared< + BuiltinFunctionEmbedder<double(const std::shared_ptr<const IDiscreteFunction>&)>>( + + [](const std::shared_ptr<const IDiscreteFunction>& c) -> double { return acoustic_dt(c); } + )); } diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index c82e9ce06ec8b7d43ce1c6cc2fec282f89dc1fee..99cc8aeeb4ca127e8e60ba3e7ee75fcc93e8db76 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -41,6 +41,7 @@ class MeshData : public IMeshData CellValue<const Rd> m_cell_centroid; CellValue<const Rd> m_cell_iso_barycenter; CellValue<const double> m_Vj; + CellValue<const double> m_sum_r_ljr; FaceValue<const double> m_ll; PUGS_INLINE @@ -435,6 +436,34 @@ class MeshData : public IMeshData } } + PUGS_INLINE + void + _computeSumOverRljr() + { + CellValue<double> sum_r_ljr(m_mesh.connectivity()); + + if constexpr (Dimension == 1) { + parallel_for( + m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { sum_r_ljr[j] = 2; }); + } else { + auto ljr = this->ljr(); + + const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix(); + + parallel_for( + m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + double sum = 0; + for (size_t r = 0; r < cell_nodes.size(); ++r) { + sum += ljr(j, r); + } + sum_r_ljr[j] = sum; + }); + } + + m_sum_r_ljr = sum_r_ljr; + } + void _checkCellVolume() { @@ -563,6 +592,16 @@ class MeshData : public IMeshData return m_Vj; } + PUGS_INLINE + CellValue<const double> + sumOverRLjr() + { + if (not m_sum_r_ljr.isBuilt()) { + this->_computeSumOverRljr(); + } + return m_sum_r_ljr; + } + private: // MeshData **must** be constructed through MeshDataManager friend class MeshDataManager; diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp new file mode 100644 index 0000000000000000000000000000000000000000..05a8adaf5f212e548e4a286ecd36f7295ad2b13b --- /dev/null +++ b/src/scheme/AcousticSolver.cpp @@ -0,0 +1,713 @@ +#include <scheme/AcousticSolver.hpp> + +#include <mesh/MeshNodeBoundary.hpp> +#include <scheme/DirichletBoundaryConditionDescriptor.hpp> +#include <scheme/DiscreteFunctionP0.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> +#include <scheme/IBoundaryConditionDescriptor.hpp> +#include <scheme/IDiscreteFunction.hpp> +#include <scheme/IDiscreteFunctionDescriptor.hpp> +#include <scheme/SymmetryBoundaryConditionDescriptor.hpp> + +#include <variant> +#include <vector> + +template <size_t Dimension> +double +acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c) +{ + const std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh = + std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(c.mesh()); + + const auto Vj = MeshDataManager::instance().getMeshData(*mesh).Vj(); + const auto Sj = MeshDataManager::instance().getMeshData(*mesh).sumOverRLjr(); + + CellValue<double> local_dt{mesh->connectivity()}; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); }); + + return min(local_dt); +} + +double +acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c) +{ + if ((c->descriptor().type() != DiscreteFunctionType::P0) or (c->dataType() != ASTNodeDataType::double_t)) { + throw NormalError("invalid discrete function type"); + } + + std::shared_ptr mesh = c->mesh(); + + switch (mesh->dimension()) { + case 1: { + return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*c)); + } + case 2: { + return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*c)); + } + case 3: { + return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*c)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +template <size_t Dimension> +class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver +{ + private: + using Rdxd = TinyMatrix<Dimension>; + using Rd = TinyVector<Dimension>; + + using MeshType = Mesh<Connectivity<Dimension>>; + using MeshDataType = MeshData<Dimension>; + + using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>; + using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>; + + class PressureBoundaryCondition; + class SymmetryBoundaryCondition; + class VelocityBoundaryCondition; + + using BoundaryCondition = + std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>; + + using BoundaryConditionList = std::vector<BoundaryCondition>; + + BoundaryConditionList m_boundary_condition_list; + + NodeValue<const Rd> m_ur; + NodeValuePerCell<const Rd> m_Fjr; + + CellValue<const double> + _getRhoC(const DiscreteScalarFunction& rho, const DiscreteScalarFunction& c) + { + Assert(rho.mesh() == c.mesh()); + + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh()); + + CellValue<double> rhoc{mesh->connectivity()}; + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; }); + + return rhoc; + } + + NodeValuePerCell<const Rdxd> + _computeGlaceAjr(const MeshType& mesh, const CellValue<const double>& rhoc) + { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + const NodeValuePerCell<const TinyVector<Dimension>> Cjr = mesh_data.Cjr(); + const NodeValuePerCell<const TinyVector<Dimension>> njr = mesh_data.njr(); + + NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()}; + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const size_t& nb_nodes = Ajr.numberOfSubValues(j); + const double& rhoc_j = rhoc[j]; + for (size_t r = 0; r < nb_nodes; ++r) { + Ajr(j, r) = tensorProduct(rhoc_j * Cjr(j, r), njr(j, r)); + } + }); + + return Ajr; + } + + NodeValue<Rdxd> + _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) + { + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + + NodeValue<Rdxd> Ar{mesh.connectivity()}; + + parallel_for( + mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { + Rdxd sum = zero; + const auto& node_to_cell = node_to_cell_matrix[r]; + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r); + + for (size_t j = 0; j < node_to_cell.size(); ++j) { + const CellId J = node_to_cell[j]; + const unsigned int R = node_local_number_in_its_cells[j]; + sum += Ajr(J, R); + } + Ar[r] = sum; + }); + + return Ar; + } + + NodeValue<Rd> + _computeBr(const Mesh<Connectivity<Dimension>>& mesh, + const NodeValuePerCell<const Rdxd>& Ajr, + const DiscreteVectorFunction& u, + const DiscreteScalarFunction& p) + { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + + NodeValue<Rd> b{mesh.connectivity()}; + + parallel_for( + mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { + const auto& node_to_cell = node_to_cell_matrix[r]; + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r); + + Rd br = zero; + for (size_t j = 0; j < node_to_cell.size(); ++j) { + const CellId J = node_to_cell[j]; + const unsigned int R = node_local_number_in_its_cells[j]; + br += Ajr(J, R) * u[J] + p[J] * Cjr(J, R); + } + + b[r] = br; + }); + + return b; + } + + BoundaryConditionList + _getBCList(const std::shared_ptr<const MeshType>& mesh, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) + { + BoundaryConditionList bc_list; + + constexpr ItemType FaceType = [] { + if constexpr (Dimension > 1) { + return ItemType::face; + } else { + return ItemType::node; + } + }(); + + for (const auto& bc_descriptor : bc_descriptor_list) { + bool is_valid_boundary_condition = true; + + switch (bc_descriptor->type()) { + case IBoundaryConditionDescriptor::Type::symmetry: { + const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor = + dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor); + for (size_t i_ref_face_list = 0; + i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { + const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); + const RefId& ref = ref_face_list.refId(); + if (ref == sym_bc_descriptor.boundaryDescriptor()) { + SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)}; + bc_list.push_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)}); + } + } + is_valid_boundary_condition = true; + break; + } + + case IBoundaryConditionDescriptor::Type::dirichlet: { + const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = + dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor); + if (dirichlet_bc_descriptor.name() == "velocity") { + for (size_t i_ref_face_list = 0; + i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { + const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); + const RefId& ref = ref_face_list.refId(); + if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) { + MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list}; + + const FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId(); + + const auto& node_list = mesh_node_boundary.nodeList(); + + Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>( + TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list); + + bc_list.push_back(VelocityBoundaryCondition{node_list, value_list}); + } + } + } else if (dirichlet_bc_descriptor.name() == "pressure") { + for (size_t i_ref_face_list = 0; + i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { + const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); + const RefId& ref = ref_face_list.refId(); + if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) { + const auto& face_list = ref_face_list.list(); + + const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId(); + + Array<const double> face_values = [&] { + if constexpr (Dimension == 1) { + return InterpolateItemValue<double( + TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh->xr(), face_list); + } else { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + return InterpolateItemValue<double( + TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list); + } + }(); + + bc_list.push_back(PressureBoundaryCondition{face_list, face_values}); + } + } + } else { + is_valid_boundary_condition = false; + } + break; + } + default: { + is_valid_boundary_condition = false; + } + } + if (not is_valid_boundary_condition) { + std::ostringstream error_msg; + error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver"; + throw NormalError(error_msg.str()); + } + } + + return bc_list; + } + + void _applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br); + void _applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br); + void _applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br); + + void + _applyBoundaryConditions(const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) + { + this->_applyPressureBC(mesh, br); + this->_applySymmetryBC(Ar, br); + this->_applyVelocityBC(Ar, br); + } + + NodeValue<const Rd> + _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) + { + NodeValue<Rd> u{mesh.connectivity()}; + parallel_for( + mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; }); + + return u; + } + + NodeValuePerCell<Rd> + _computeFjr(const MeshType& mesh, + const NodeValuePerCell<const Rdxd>& Ajr, + const NodeValue<const Rd>& ur, + const DiscreteVectorFunction& u, + const DiscreteScalarFunction& p) + { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + + const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); + + NodeValuePerCell<Rd> F{mesh.connectivity()}; + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_nodes = cell_to_node_matrix[j]; + + for (size_t r = 0; r < cell_nodes.size(); ++r) { + F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r); + } + }); + + return F; + } + + AcousticSolver(const std::shared_ptr<const MeshType>& p_mesh, + const DiscreteScalarFunction& rho, + const DiscreteScalarFunction& c, + const DiscreteVectorFunction& u, + const DiscreteScalarFunction& p, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) + : m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)} + { + const MeshType& mesh = *p_mesh; + + CellValue<const double> rhoc = this->_getRhoC(rho, c); + NodeValuePerCell<const Rdxd> Ajr = this->_computeGlaceAjr(mesh, rhoc); + + NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr); + NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p); + + this->_applyBoundaryConditions(mesh, Ar, br); + + synchronize(Ar); + synchronize(br); + + m_ur = this->_computeUr(mesh, Ar, br); + m_Fjr = this->_computeFjr(mesh, Ajr, m_ur, u, p); + } + + public: + std::tuple<std::shared_ptr<const IMesh>, + std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>, + std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>, + std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>> + apply(const double& dt, + const std::shared_ptr<const MeshType>& mesh, + const DiscreteFunctionP0<Dimension, double>& rho, + const DiscreteFunctionP0<Dimension, Rd>& u, + const DiscreteFunctionP0<Dimension, double>& E) const + { + const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix(); + + 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 += (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 {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<std::shared_ptr<const IMesh>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>> + apply(const double& dt, + 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 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(dt, 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)); + } + + AcousticSolver(const std::shared_ptr<const IMesh>& mesh, + const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& c, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) + : AcousticSolver{std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh), + *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho), + *std::dynamic_pointer_cast<const DiscreteScalarFunction>(c), + *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u), + *std::dynamic_pointer_cast<const DiscreteScalarFunction>(p), + bc_descriptor_list} + {} + + AcousticSolver() = default; + AcousticSolver(AcousticSolver&&) = default; + ~AcousticSolver() = default; +}; + +template <size_t Dimension> +void +AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br) +{ + for (const auto& boundary_condition : m_boundary_condition_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<PressureBoundaryCondition, T>) { + MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + if constexpr (Dimension == 1) { + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + + const auto& node_list = bc.faceList(); + const auto& value_list = bc.valueList(); + parallel_for( + node_list.size(), PUGS_LAMBDA(size_t i_node) { + const NodeId node_id = node_list[i_node]; + const auto& node_cell_list = node_to_cell_matrix[node_id]; + Assert(node_cell_list.size() == 1); + + CellId node_cell_id = node_cell_list[0]; + size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0); + + br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell); + }); + } else { + const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr(); + + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& value_list = bc.valueList(); + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + const auto& face_cell_list = face_to_cell_matrix[face_id]; + Assert(face_cell_list.size() == 1); + + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0); + + const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1; + + const auto& face_nodes = face_to_node_matrix[face_id]; + + for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) { + NodeId node_id = face_nodes[i_node]; + br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node); + } + } + } + } + }, + boundary_condition); + } +} + +template <size_t Dimension> +void +AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) +{ + for (const auto& boundary_condition : m_boundary_condition_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) { + const Rd& n = bc.outgoingNormal(); + + const Rdxd I = identity; + const Rdxd nxn = tensorProduct(n, n); + const Rdxd P = I - nxn; + + const Array<const NodeId>& node_list = bc.nodeList(); + parallel_for( + bc.numberOfNodes(), PUGS_LAMBDA(int r_number) { + const NodeId r = node_list[r_number]; + + Ar[r] = P * Ar[r] * P + nxn; + br[r] = P * br[r]; + }); + } + }, + boundary_condition); + } +} + +template <size_t Dimension> +void +AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) +{ + for (const auto& boundary_condition : m_boundary_condition_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) { + const auto& node_list = bc.nodeList(); + const auto& value_list = bc.valueList(); + + parallel_for( + node_list.size(), PUGS_LAMBDA(size_t i_node) { + NodeId node_id = node_list[i_node]; + const auto& value = value_list[i_node]; + + Ar[node_id] = identity; + br[node_id] = value; + }); + } + }, + boundary_condition); + } +} + +template <size_t Dimension> +class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryCondition +{ + private: + const Array<const double> m_value_list; + const Array<const FaceId> m_face_list; + + public: + const Array<const FaceId>& + faceList() const + { + return m_face_list; + } + + const Array<const double>& + valueList() const + { + return m_value_list; + } + + PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list) + : m_value_list{value_list}, m_face_list{face_list} + {} + + ~PressureBoundaryCondition() = default; +}; + +template <> +class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition +{ + private: + const Array<const double> m_value_list; + const Array<const NodeId> m_face_list; + + public: + const Array<const NodeId>& + faceList() const + { + return m_face_list; + } + + const Array<const double>& + valueList() const + { + return m_value_list; + } + + PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list) + : m_value_list{value_list}, m_face_list{face_list} + {} + + ~PressureBoundaryCondition() = default; +}; + +template <size_t Dimension> +class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryCondition +{ + private: + const Array<const TinyVector<Dimension>> m_value_list; + const Array<const NodeId> m_node_list; + + public: + const Array<const NodeId>& + nodeList() const + { + return m_node_list; + } + + const Array<const TinyVector<Dimension>>& + valueList() const + { + return m_value_list; + } + + VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list) + : m_value_list{value_list}, m_node_list{node_list} + {} + + ~VelocityBoundaryCondition() = default; +}; + +template <size_t Dimension> +class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary) + { + ; + } + + ~SymmetryBoundaryCondition() = default; +}; + +AcousticSolverHandler::AcousticSolverHandler( + const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& c, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) +{ + std::shared_ptr i_mesh = getCommonMesh({rho, c, u, p}); + if (not i_mesh) { + throw NormalError("discrete functions are not defined on the same mesh"); + } + + if (not checkDiscretizationType({rho, c, u, p}, DiscreteFunctionType::P0)) { + throw NormalError("acoustic solver expects P0 functions"); + } + + switch (i_mesh->dimension()) { + case 1: { + m_acoustic_solver = std::make_unique<AcousticSolver<1>>(i_mesh, rho, c, u, p, bc_descriptor_list); + break; + } + case 2: { + m_acoustic_solver = std::make_unique<AcousticSolver<2>>(i_mesh, rho, c, u, p, bc_descriptor_list); + break; + } + case 3: { + m_acoustic_solver = std::make_unique<AcousticSolver<3>>(i_mesh, rho, c, u, p, bc_descriptor_list); + break; + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 3e3151d25429593955cb281c2dd814690f67bfa5..6618361cc56a6efe55620fda83c8976ecd60b3f1 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -21,8 +21,53 @@ #include <iostream> +class IDiscreteFunction; +class IBoundaryConditionDescriptor; + +double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c); + +class AcousticSolverHandler +{ + private: + struct IAcousticSolver + { + virtual std::tuple<std::shared_ptr<const IMesh>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>, + std::shared_ptr<const IDiscreteFunction>> + apply(const double& dt, + const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& E) const = 0; + + IAcousticSolver() = default; + IAcousticSolver(IAcousticSolver&&) = default; + IAcousticSolver& operator=(IAcousticSolver&&) = default; + + virtual ~IAcousticSolver() = default; + }; + + template <size_t Dimension> + class AcousticSolver; + + std::unique_ptr<IAcousticSolver> m_acoustic_solver; + + public: + const IAcousticSolver& + solver() const + { + return *m_acoustic_solver; + } + + AcousticSolverHandler(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& c, + const std::shared_ptr<const IDiscreteFunction>& u, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list); +}; + template <typename MeshType> -class AcousticSolver +class AcousticSolverOld { public: class PressureBoundaryCondition; @@ -392,9 +437,9 @@ class AcousticSolver CellValue<double> m_Vj_over_cj; public: - AcousticSolver(std::shared_ptr<const MeshType> p_mesh, - const BoundaryConditionList& bc_list, - const AcousticSolverType solver_type) + AcousticSolverOld(std::shared_ptr<const MeshType> p_mesh, + const BoundaryConditionList& bc_list, + const AcousticSolverType solver_type) : m_mesh(p_mesh), m_connectivity(m_mesh->connectivity()), m_boundary_condition_list(bc_list), @@ -490,7 +535,7 @@ class AcousticSolver }; template <typename MeshType> -class AcousticSolver<MeshType>::PressureBoundaryCondition +class AcousticSolverOld<MeshType>::PressureBoundaryCondition { private: const Array<const double> m_value_list; @@ -517,7 +562,7 @@ class AcousticSolver<MeshType>::PressureBoundaryCondition }; template <> -class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition +class AcousticSolverOld<Mesh<Connectivity<1>>>::PressureBoundaryCondition { private: const Array<const double> m_value_list; @@ -544,7 +589,7 @@ class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition }; template <typename MeshType> -class AcousticSolver<MeshType>::VelocityBoundaryCondition +class AcousticSolverOld<MeshType>::VelocityBoundaryCondition { private: const Array<const TinyVector<Dimension>> m_value_list; @@ -571,7 +616,7 @@ class AcousticSolver<MeshType>::VelocityBoundaryCondition }; template <typename MeshType> -class AcousticSolver<MeshType>::SymmetryBoundaryCondition +class AcousticSolverOld<MeshType>::SymmetryBoundaryCondition { public: static constexpr size_t Dimension = MeshType::Dimension; diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt index 1467393f2cbcbd9e31dd495e2905e9be478a1d78..b7522b9bbd08a90ee1bd95a56747753d902e10ed 100644 --- a/src/scheme/CMakeLists.txt +++ b/src/scheme/CMakeLists.txt @@ -2,4 +2,7 @@ add_library( PugsScheme - DiscreteFunctionInterpoler.cpp) + AcousticSolver.cpp + DiscreteFunctionInterpoler.cpp + DiscreteFunctionUtils.cpp + PerfectGas.cpp) diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp index bf97079b776c80fec10e26608cfc6ff5e40a664a..c51e8e09dbb78dc6f40f4fbb18eb737ddb7e1a30 100644 --- a/src/scheme/DiscreteFunctionP0.hpp +++ b/src/scheme/DiscreteFunctionP0.hpp @@ -28,7 +28,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction return ast_node_data_type_from<DataType>; } - CellValue<DataType> + const CellValue<DataType>& cellValues() const { return m_cell_values; @@ -63,6 +63,17 @@ class DiscreteFunctionP0 : public IDiscreteFunction mesh_data.xj()); } + DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()} + { + ; + } + + DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value) + : m_mesh{mesh}, m_cell_values{cell_value} + { + ; + } + DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default; DiscreteFunctionP0(DiscreteFunctionP0&&) noexcept = default; diff --git a/src/scheme/DiscreteFunctionUtils.cpp b/src/scheme/DiscreteFunctionUtils.cpp new file mode 100644 index 0000000000000000000000000000000000000000..56453a746a19b9ee999fe1548d7067250e2884eb --- /dev/null +++ b/src/scheme/DiscreteFunctionUtils.cpp @@ -0,0 +1,115 @@ +#include <scheme/DiscreteFunctionUtils.hpp> + +#include <mesh/Connectivity.hpp> +#include <mesh/IMesh.hpp> +#include <mesh/Mesh.hpp> +#include <scheme/DiscreteFunctionP0.hpp> + +template <size_t Dimension, typename DataType> +std::shared_ptr<const IDiscreteFunction> +shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh, + const std::shared_ptr<const DiscreteFunctionP0<Dimension, DataType>>& discrete_function) +{ + return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh, discrete_function->cellValues()); +} + +template <size_t Dimension> +std::shared_ptr<const IDiscreteFunction> +shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh, + const std::shared_ptr<const IDiscreteFunction>& discrete_function) +{ + const std::shared_ptr function_mesh = + std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(discrete_function->mesh()); + + if (mesh->shared_connectivity() != function_mesh->shared_connectivity()) { + throw NormalError("incompatible connectivities"); + } + + switch (discrete_function->descriptor().type()) { + case DiscreteFunctionType::P0: { + switch (discrete_function->dataType()) { + case ASTNodeDataType::double_t: { + return shallowCopy(mesh, + std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, double>>(discrete_function)); + } + case ASTNodeDataType::vector_t: { + switch (discrete_function->dataType().dimension()) { + case 1: { + return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>>( + discrete_function)); + } + case 2: { + return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<2>>>( + discrete_function)); + } + case 3: { + return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>>( + discrete_function)); + } + default: { + throw UnexpectedError("invalid data vector dimension: " + + std::to_string(discrete_function->dataType().dimension())); + } + } + } + case ASTNodeDataType::matrix_t: { + if (discrete_function->dataType().nbRows() != discrete_function->dataType().nbColumns()) { + throw UnexpectedError( + "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" + + std::to_string(discrete_function->dataType().nbColumns())); + } + switch (discrete_function->dataType().nbRows()) { + case 1: { + return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>( + discrete_function)); + } + case 2: { + return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>( + discrete_function)); + } + case 3: { + return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>( + discrete_function)); + } + default: { + throw UnexpectedError( + "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" + + std::to_string(discrete_function->dataType().nbColumns())); + } + } + } + default: { + throw UnexpectedError("invalid kind of P0 function: invalid data type"); + } + } + } + default: { + throw NormalError("invalid discretization type"); + } + } +} + +std::shared_ptr<const IDiscreteFunction> +shallowCopy(const std::shared_ptr<const IMesh>& mesh, const std::shared_ptr<const IDiscreteFunction>& discrete_function) +{ + if (mesh == discrete_function->mesh()) { + return discrete_function; + } else if (mesh->dimension() != discrete_function->mesh()->dimension()) { + throw NormalError("incompatible mesh dimensions"); + } + + switch (mesh->dimension()) { + case 1: { + return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), discrete_function); + } + case 2: { + return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), discrete_function); + } + case 3: { + return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), discrete_function); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} diff --git a/src/scheme/DiscreteFunctionUtils.hpp b/src/scheme/DiscreteFunctionUtils.hpp new file mode 100644 index 0000000000000000000000000000000000000000..1eea93a7eae8db4c592b6bbec35341ad33e35bdc --- /dev/null +++ b/src/scheme/DiscreteFunctionUtils.hpp @@ -0,0 +1,43 @@ +#ifndef DISCRETE_FUNCTION_UTILS_HPP +#define DISCRETE_FUNCTION_UTILS_HPP + +#include <scheme/DiscreteFunctionType.hpp> +#include <scheme/IDiscreteFunction.hpp> +#include <scheme/IDiscreteFunctionDescriptor.hpp> + +#include <vector> + +PUGS_INLINE +bool +checkDiscretizationType(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list, + const DiscreteFunctionType& discrete_function_type) +{ + for (auto discrete_function : discrete_function_list) { + if (discrete_function->descriptor().type() != discrete_function_type) { + return false; + } + } + return true; +} + +PUGS_INLINE +std::shared_ptr<const IMesh> +getCommonMesh(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list) +{ + std::shared_ptr<const IMesh> i_mesh; + for (auto discrete_function : discrete_function_list) { + if (not i_mesh.use_count()) { + i_mesh = discrete_function->mesh(); + } else { + if (i_mesh != discrete_function->mesh()) { + return {}; + } + } + } + return i_mesh; +} + +std::shared_ptr<const IDiscreteFunction> shallowCopy(const std::shared_ptr<const IMesh>& mesh, + const std::shared_ptr<const IDiscreteFunction>& discrete_function); + +#endif // DISCRETE_FUNCTION_UTILS_HPP diff --git a/src/scheme/PerfectGas.cpp b/src/scheme/PerfectGas.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b371dea3c9e4aac8ee867dd0f938b72bcd02791e --- /dev/null +++ b/src/scheme/PerfectGas.cpp @@ -0,0 +1,248 @@ +#include <scheme/PerfectGas.hpp> + +#include <scheme/DiscreteFunctionP0.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> + +namespace perfect_gas +{ +template <size_t Dimension, typename DataType> +std::shared_ptr<IDiscreteFunction> +soundSpeed(const DiscreteFunctionP0<Dimension, DataType>& rho, + const DiscreteFunctionP0<Dimension, DataType>& p, + const DiscreteFunctionP0<Dimension, DataType>& gamma) +{ + using MeshType = Mesh<Connectivity<Dimension>>; + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh()); + + DiscreteFunctionP0<Dimension, DataType> sound_speed{mesh}; + + parallel_for( + mesh->numberOfCells(), + PUGS_LAMBDA(CellId cell_id) { sound_speed[cell_id] = std::sqrt(gamma[cell_id] * p[cell_id] / rho[cell_id]); }); + + return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(sound_speed); +} + +std::shared_ptr<IDiscreteFunction> +soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::shared_ptr<const IDiscreteFunction>& gamma) +{ + std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma}); + if (not mesh.use_count()) { + throw NormalError("rho, p and gamma are not defined on the same mesh"); + } + + switch (mesh->dimension()) { + case 1: { + return soundSpeed(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p), + dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma)); + } + case 2: { + return soundSpeed(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p), + dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma)); + } + case 3: { + return soundSpeed(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p), + dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +template <size_t Dimension, typename DataType> +std::shared_ptr<IDiscreteFunction> +epsilonFromRhoPGamma(const DiscreteFunctionP0<Dimension, DataType>& rho, + const DiscreteFunctionP0<Dimension, DataType>& p, + const DiscreteFunctionP0<Dimension, DataType>& gamma) +{ + using MeshType = Mesh<Connectivity<Dimension>>; + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh()); + + DiscreteFunctionP0<Dimension, DataType> epsilon{mesh}; + + parallel_for( + mesh->numberOfCells(), + PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = p[cell_id] / ((gamma[cell_id] - 1) * rho[cell_id]); }); + + return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon); +} + +std::shared_ptr<IDiscreteFunction> +epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::shared_ptr<const IDiscreteFunction>& gamma) +{ + std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma}); + if (not mesh.use_count()) { + throw NormalError("rho, p and gamma are not defined on the same mesh"); + } + + switch (mesh->dimension()) { + case 1: { + return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p), + dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma)); + } + case 2: { + return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p), + dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma)); + } + case 3: { + return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p), + dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +template <size_t Dimension, typename DataType> +std::shared_ptr<IDiscreteFunction> +pFromRhoEpsilonGamma(const DiscreteFunctionP0<Dimension, DataType>& rho, + const DiscreteFunctionP0<Dimension, DataType>& epsilon, + const DiscreteFunctionP0<Dimension, DataType>& gamma) +{ + using MeshType = Mesh<Connectivity<Dimension>>; + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh()); + + DiscreteFunctionP0<Dimension, DataType> p{mesh}; + + parallel_for( + mesh->numberOfCells(), + PUGS_LAMBDA(CellId cell_id) { p[cell_id] = ((gamma[cell_id] - 1) * rho[cell_id]) * epsilon[cell_id]; }); + + return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(p); +} + +std::shared_ptr<IDiscreteFunction> +pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& epsilon, + const std::shared_ptr<const IDiscreteFunction>& gamma) +{ + std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, epsilon, gamma}); + if (not mesh.use_count()) { + throw NormalError("rho, epsilon and gamma are not defined on the same mesh"); + } + + switch (mesh->dimension()) { + case 1: { + return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon), + dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma)); + } + case 2: { + return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon), + dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma)); + } + case 3: { + return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho), + dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon), + dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +template <size_t Dimension, typename DataType> +std::shared_ptr<IDiscreteFunction> +totalEnergyFromEpsilonU(const DiscreteFunctionP0<Dimension, DataType>& epsilon, + const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u) +{ + using MeshType = Mesh<Connectivity<Dimension>>; + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(epsilon.mesh()); + + DiscreteFunctionP0<Dimension, DataType> E{mesh}; + + parallel_for( + mesh->numberOfCells(), + PUGS_LAMBDA(CellId cell_id) { E[cell_id] = epsilon[cell_id] + 0.5 * (u[cell_id], u[cell_id]); }); + + return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon); +} + +std::shared_ptr<IDiscreteFunction> +totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon, + const std::shared_ptr<const IDiscreteFunction>& u) +{ + std::shared_ptr<const IMesh> mesh = getCommonMesh({epsilon, u}); + if (not mesh.use_count()) { + throw NormalError("epsilon and u are not defined on the same mesh"); + } + + switch (mesh->dimension()) { + case 1: { + return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon), + dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u)); + } + case 2: { + return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon), + dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u)); + } + case 3: { + return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon), + dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +template <size_t Dimension, typename DataType> +std::shared_ptr<IDiscreteFunction> +epsilonFromTotalEnergyU(const DiscreteFunctionP0<Dimension, DataType>& E, + const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u) +{ + using MeshType = Mesh<Connectivity<Dimension>>; + std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(u.mesh()); + + DiscreteFunctionP0<Dimension, DataType> epsilon{mesh}; + + parallel_for( + mesh->numberOfCells(), + PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = E[cell_id] - 0.5 * (u[cell_id], u[cell_id]); }); + + return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon); +} + +std::shared_ptr<IDiscreteFunction> +epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E, + const std::shared_ptr<const IDiscreteFunction>& u) +{ + std::shared_ptr<const IMesh> mesh = getCommonMesh({E, u}); + if (not mesh.use_count()) { + throw NormalError("E and u are not defined on the same mesh"); + } + + switch (mesh->dimension()) { + case 1: { + return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*E), + dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u)); + } + case 2: { + return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*E), + dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u)); + } + case 3: { + return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*E), + dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u)); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +} // namespace perfect_gas diff --git a/src/scheme/PerfectGas.hpp b/src/scheme/PerfectGas.hpp new file mode 100644 index 0000000000000000000000000000000000000000..fb7b9d212b4681f365faa04e4565dfd78c316f66 --- /dev/null +++ b/src/scheme/PerfectGas.hpp @@ -0,0 +1,32 @@ +#ifndef PERFECT_GAS_HPP +#define PERFECT_GAS_HPP + +#include <scheme/IDiscreteFunction.hpp> + +namespace perfect_gas +{ +std::shared_ptr<IDiscreteFunction> soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::shared_ptr<const IDiscreteFunction>& gamma); + +std::shared_ptr<IDiscreteFunction> epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& p, + const std::shared_ptr<const IDiscreteFunction>& gamma); + +std::shared_ptr<IDiscreteFunction> pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho, + const std::shared_ptr<const IDiscreteFunction>& epsilon, + const std::shared_ptr<const IDiscreteFunction>& gamma); + +// This is a temporary function +// Todo: TO REMOVE +std::shared_ptr<IDiscreteFunction> totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon, + const std::shared_ptr<const IDiscreteFunction>& u); + +// This is a temporary function +// Todo: TO REMOVE +std::shared_ptr<IDiscreteFunction> epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E, + const std::shared_ptr<const IDiscreteFunction>& u); + +} // namespace perfect_gas + +#endif // PERFECT_GAS_HPP