From a3c94db15b186fb4513b1483c78bf3d6e3789342 Mon Sep 17 00:00:00 2001 From: Stephane Del Pino <stephane.delpino44@gmail.com> Date: Tue, 28 Jul 2020 19:02:43 +0200 Subject: [PATCH] Prepare compilation time reduction --- .../algorithms/AcousticSolverAlgorithm.hpp | 222 +++++++++++++ src/language/modules/SchemeModule.cpp | 304 +----------------- src/language/utils/InterpolateItemValue.hpp | 82 +++++ 3 files changed, 311 insertions(+), 297 deletions(-) create mode 100644 src/language/algorithms/AcousticSolverAlgorithm.hpp create mode 100644 src/language/utils/InterpolateItemValue.hpp diff --git a/src/language/algorithms/AcousticSolverAlgorithm.hpp b/src/language/algorithms/AcousticSolverAlgorithm.hpp new file mode 100644 index 000000000..3fca4703f --- /dev/null +++ b/src/language/algorithms/AcousticSolverAlgorithm.hpp @@ -0,0 +1,222 @@ +#ifndef ACOUSTIC_SOLVER_ALGORITHM_HPP +#define ACOUSTIC_SOLVER_ALGORITHM_HPP + +#include <language/utils/FunctionSymbolId.hpp> +#include <language/utils/InterpolateItemValue.hpp> +#include <output/VTKWriter.hpp> +#include <scheme/AcousticSolver.hpp> +#include <scheme/IBoundaryConditionDescriptor.hpp> +#include <scheme/IBoundaryDescriptor.hpp> +#include <scheme/PressureBoundaryConditionDescriptor.hpp> +#include <scheme/SymmetryBoundaryConditionDescriptor.hpp> +#include <scheme/VelocityBoundaryConditionDescriptor.hpp> + +template <size_t Dimension, AcousticSolverType solver_type> +struct AcousticSolverAlgorithm +{ + using ConnectivityType = Connectivity<Dimension>; + using MeshType = Mesh<ConnectivityType>; + using MeshDataType = MeshData<Dimension>; + using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>; + + AcousticSolverAlgorithm(std::shared_ptr<const IMesh> i_mesh, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const FunctionSymbolId& rho_id, + const FunctionSymbolId& u_id, + const FunctionSymbolId& p_id) + { + std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh); + + std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n'; + + typename AcousticSolver<MeshType>::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) { + switch (bc_descriptor->type()) { + case IBoundaryConditionDescriptor::Type::symmetry: { + using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition; + + 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()) { + bc_list.push_back( + SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list)}); + } + } + break; + } + case IBoundaryConditionDescriptor::Type::velocity: { + using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition; + + const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor = + dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*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 == velocity_bc_descriptor.boundaryDescriptor()) { + MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list}; + + const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId(); + + 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}); + } + } + break; + } + case IBoundaryConditionDescriptor::Type::pressure: { + using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition; + + const PressureBoundaryConditionDescriptor& pressure_bc_descriptor = + dynamic_cast<const PressureBoundaryConditionDescriptor&>(*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 == pressure_bc_descriptor.boundaryDescriptor()) { + const auto& face_list = ref_face_list.list(); + + const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId(); + + 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}); + } + } + break; + } + default: { + std::ostringstream error_msg; + error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver"; + throw NormalError(error_msg.str()); + } + } + } + } + + UnknownsType unknowns(*mesh); + + { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + unknowns.rhoj() = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, + mesh_data.xj()); + + unknowns.pj() = + InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj()); + + unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>( + TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj()); + } + unknowns.gammaj().fill(1.4); + + AcousticSolver acoustic_solver(mesh, bc_list, solver_type); + + const double tmax = 0.2; + double t = 0; + + int itermax = std::numeric_limits<int>::max(); + int iteration = 0; + + CellValue<double>& rhoj = unknowns.rhoj(); + CellValue<double>& ej = unknowns.ej(); + CellValue<double>& pj = unknowns.pj(); + CellValue<double>& gammaj = unknowns.gammaj(); + CellValue<double>& cj = unknowns.cj(); + CellValue<TinyVector<Dimension>>& uj = unknowns.uj(); + CellValue<double>& Ej = unknowns.Ej(); + CellValue<double>& mj = unknowns.mj(); + CellValue<double>& inv_mj = unknowns.invMj(); + + BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); + block_eos.updateEandCFromRhoP(); + + { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + const CellValue<const double> Vj = mesh_data.Vj(); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); }); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; }); + + parallel_for( + mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; }); + } + + VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01); + + while ((t < tmax) and (iteration < itermax)) { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + vtk_writer.write(mesh, + {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()}, + NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, + NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()}, + NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}}, + t); + + double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj); + if (t + dt > tmax) { + dt = tmax - t; + } + + std::cout.setf(std::cout.scientific); + std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset + << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt + << rang::style::reset << '\n'; + + mesh = acoustic_solver.computeNextStep(dt, unknowns); + + block_eos.updatePandCFromRhoE(); + + t += dt; + ++iteration; + } + std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after " + << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations" + << rang::style::reset << '\n'; + { + MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh); + + vtk_writer.write(mesh, + {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()}, + NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, + NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()}, + NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}}, + t, true); // forces last output + } + } +}; + +#endif // ACOUSTIC_SOLVER_ALGORITHM_HPP diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 477e3cc80..b38373e67 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -1,5 +1,6 @@ #include <language/modules/SchemeModule.hpp> +#include <language/algorithms/AcousticSolverAlgorithm.hpp> #include <language/utils/BuiltinFunctionEmbedder.hpp> #include <language/utils/TypeDescriptor.hpp> #include <mesh/Mesh.hpp> @@ -8,300 +9,9 @@ #include <scheme/IBoundaryDescriptor.hpp> #include <scheme/NamedBoundaryDescriptor.hpp> #include <scheme/NumberedBoundaryDescriptor.hpp> -#include <scheme/PressureBoundaryConditionDescriptor.hpp> -#include <scheme/SymmetryBoundaryConditionDescriptor.hpp> -#include <scheme/VelocityBoundaryConditionDescriptor.hpp> #include <memory> -/////////// TEMPORARY - -#include <language/utils/PugsFunctionAdapter.hpp> -#include <output/VTKWriter.hpp> - -template <typename T> -class InterpolateItemValue; -template <typename OutputType, typename InputType> -class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)> -{ - static constexpr size_t Dimension = OutputType::Dimension; - using Adapter = PugsFunctionAdapter<OutputType(InputType)>; - - public: - template <ItemType item_type> - static inline ItemValue<OutputType, item_type> - interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position) - { - auto& expression = Adapter::getFunctionExpression(function_symbol_id); - auto convert_result = Adapter::getResultConverter(expression.m_data_type); - - Array<ExecutionPolicy> context_list = Adapter::getContextList(expression); - - using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space; - Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens; - const IConnectivity& connectivity = *position.connectivity_ptr(); - - ItemValue<OutputType, item_type> value(connectivity); - using ItemId = ItemIdT<item_type>; - - parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) { - const int32_t t = tokens.acquire(); - - auto& execution_policy = context_list[t]; - - Adapter::convertArgs(execution_policy.currentContext(), position[i]); - auto result = expression.execute(execution_policy); - value[i] = convert_result(std::move(result)); - - tokens.release(t); - }); - - return value; - } - - template <ItemType item_type> - static inline Array<OutputType> - interpolate(const FunctionSymbolId& function_symbol_id, - const ItemValue<const InputType, item_type>& position, - const Array<const ItemIdT<item_type>>& list_of_items) - { - auto& expression = Adapter::getFunctionExpression(function_symbol_id); - auto convert_result = Adapter::getResultConverter(expression.m_data_type); - - Array<ExecutionPolicy> context_list = Adapter::getContextList(expression); - - using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space; - Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens; - - Array<OutputType> value{list_of_items.size()}; - using ItemId = ItemIdT<item_type>; - - parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) { - ItemId item_id = list_of_items[i_item]; - const int32_t t = tokens.acquire(); - - auto& execution_policy = context_list[t]; - - Adapter::convertArgs(execution_policy.currentContext(), position[item_id]); - auto result = expression.execute(execution_policy); - value[i_item] = convert_result(std::move(result)); - - tokens.release(t); - }); - - return value; - } -}; - -template <size_t Dimension, AcousticSolverType solver_type> -struct AcousticSolverScheme -{ - using ConnectivityType = Connectivity<Dimension>; - using MeshType = Mesh<ConnectivityType>; - using MeshDataType = MeshData<Dimension>; - using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>; - - std::shared_ptr<const MeshType> m_mesh; - - AcousticSolverScheme(std::shared_ptr<const IMesh> i_mesh, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, - const FunctionSymbolId& rho_id, - const FunctionSymbolId& u_id, - const FunctionSymbolId& p_id) - : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)} - { - std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n'; - - typename AcousticSolver<MeshType>::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) { - switch (bc_descriptor->type()) { - case IBoundaryConditionDescriptor::Type::symmetry: { - using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition; - - const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor = - dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor); - for (size_t i_ref_face_list = 0; - i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { - const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); - const RefId& ref = ref_face_list.refId(); - if (ref == sym_bc_descriptor.boundaryDescriptor()) { - bc_list.push_back( - SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list)}); - } - } - break; - } - case IBoundaryConditionDescriptor::Type::velocity: { - using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition; - - const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor = - dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*bc_descriptor); - for (size_t i_ref_face_list = 0; - i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { - const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); - const RefId& ref = ref_face_list.refId(); - if (ref == velocity_bc_descriptor.boundaryDescriptor()) { - MeshNodeBoundary<Dimension> mesh_node_boundary{m_mesh, ref_face_list}; - - const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId(); - - 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, m_mesh->xr(), node_list); - - bc_list.push_back(VelocityBoundaryCondition{node_list, value_list}); - } - } - break; - } - case IBoundaryConditionDescriptor::Type::pressure: { - using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition; - - const PressureBoundaryConditionDescriptor& pressure_bc_descriptor = - dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor); - for (size_t i_ref_face_list = 0; - i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) { - const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list); - const RefId& ref = ref_face_list.refId(); - if (ref == pressure_bc_descriptor.boundaryDescriptor()) { - const auto& face_list = ref_face_list.list(); - - const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId(); - - Array<const double> face_values = [&] { - if constexpr (Dimension == 1) { - return InterpolateItemValue<double( - TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, m_mesh->xr(), face_list); - } else { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_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}); - } - } - break; - } - default: { - std::ostringstream error_msg; - error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver"; - throw NormalError(error_msg.str()); - } - } - } - } - - UnknownsType unknowns(*m_mesh); - - { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); - - unknowns.rhoj() = - InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, - mesh_data.xj()); - - unknowns.pj() = - InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj()); - - unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>( - TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj()); - } - unknowns.gammaj().fill(1.4); - - AcousticSolver acoustic_solver(m_mesh, bc_list, solver_type); - - const double tmax = 0.2; - double t = 0; - - int itermax = std::numeric_limits<int>::max(); - int iteration = 0; - - CellValue<double>& rhoj = unknowns.rhoj(); - CellValue<double>& ej = unknowns.ej(); - CellValue<double>& pj = unknowns.pj(); - CellValue<double>& gammaj = unknowns.gammaj(); - CellValue<double>& cj = unknowns.cj(); - CellValue<TinyVector<Dimension>>& uj = unknowns.uj(); - CellValue<double>& Ej = unknowns.Ej(); - CellValue<double>& mj = unknowns.mj(); - CellValue<double>& inv_mj = unknowns.invMj(); - - BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - block_eos.updateEandCFromRhoP(); - - { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); - - const CellValue<const double> Vj = mesh_data.Vj(); - - parallel_for( - m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); }); - - parallel_for( - m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; }); - - parallel_for( - m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; }); - } - - VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01); - - while ((t < tmax) and (iteration < itermax)) { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); - - vtk_writer.write(m_mesh, - {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()}, - NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, - NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()}, - NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}}, - t); - - double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj); - if (t + dt > tmax) { - dt = tmax - t; - } - - std::cout.setf(std::cout.scientific); - std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset - << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt - << rang::style::reset << '\n'; - - m_mesh = acoustic_solver.computeNextStep(dt, unknowns); - - block_eos.updatePandCFromRhoE(); - - t += dt; - ++iteration; - } - std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after " - << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations" - << rang::style::reset << '\n'; - { - MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh); - - vtk_writer.write(m_mesh, - {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()}, - NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()}, - NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()}, - NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}}, - t, true); // forces last output - } - } -}; - SchemeModule::SchemeModule() { this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>); @@ -376,15 +86,15 @@ SchemeModule::SchemeModule() const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void { switch (p_mesh->dimension()) { case 1: { - AcousticSolverScheme<1, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; + AcousticSolverAlgorithm<1, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; break; } case 2: { - AcousticSolverScheme<2, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; + AcousticSolverAlgorithm<2, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; break; } case 3: { - AcousticSolverScheme<3, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; + AcousticSolverAlgorithm<3, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; break; } default: { @@ -406,15 +116,15 @@ SchemeModule::SchemeModule() const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void { switch (p_mesh->dimension()) { case 1: { - AcousticSolverScheme<1, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; + AcousticSolverAlgorithm<1, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; break; } case 2: { - AcousticSolverScheme<2, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; + AcousticSolverAlgorithm<2, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; break; } case 3: { - AcousticSolverScheme<3, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; + AcousticSolverAlgorithm<3, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id}; break; } default: { diff --git a/src/language/utils/InterpolateItemValue.hpp b/src/language/utils/InterpolateItemValue.hpp new file mode 100644 index 000000000..16a87b495 --- /dev/null +++ b/src/language/utils/InterpolateItemValue.hpp @@ -0,0 +1,82 @@ +#ifndef INTERPOLATE_ITEM_VALUE_HPP +#define INTERPOLATE_ITEM_VALUE_HPP + +#include <language/utils/PugsFunctionAdapter.hpp> +#include <mesh/ItemType.hpp> +#include <mesh/ItemValue.hpp> + +template <typename T> +class InterpolateItemValue; +template <typename OutputType, typename InputType> +class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)> +{ + static constexpr size_t Dimension = OutputType::Dimension; + using Adapter = PugsFunctionAdapter<OutputType(InputType)>; + + public: + template <ItemType item_type> + static inline ItemValue<OutputType, item_type> + interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position) + { + auto& expression = Adapter::getFunctionExpression(function_symbol_id); + auto convert_result = Adapter::getResultConverter(expression.m_data_type); + + Array<ExecutionPolicy> context_list = Adapter::getContextList(expression); + + using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space; + Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens; + const IConnectivity& connectivity = *position.connectivity_ptr(); + + ItemValue<OutputType, item_type> value(connectivity); + using ItemId = ItemIdT<item_type>; + + parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) { + const int32_t t = tokens.acquire(); + + auto& execution_policy = context_list[t]; + + Adapter::convertArgs(execution_policy.currentContext(), position[i]); + auto result = expression.execute(execution_policy); + value[i] = convert_result(std::move(result)); + + tokens.release(t); + }); + + return value; + } + + template <ItemType item_type> + static inline Array<OutputType> + interpolate(const FunctionSymbolId& function_symbol_id, + const ItemValue<const InputType, item_type>& position, + const Array<const ItemIdT<item_type>>& list_of_items) + { + auto& expression = Adapter::getFunctionExpression(function_symbol_id); + auto convert_result = Adapter::getResultConverter(expression.m_data_type); + + Array<ExecutionPolicy> context_list = Adapter::getContextList(expression); + + using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space; + Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens; + + Array<OutputType> value{list_of_items.size()}; + using ItemId = ItemIdT<item_type>; + + parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) { + ItemId item_id = list_of_items[i_item]; + const int32_t t = tokens.acquire(); + + auto& execution_policy = context_list[t]; + + Adapter::convertArgs(execution_policy.currentContext(), position[item_id]); + auto result = expression.execute(execution_policy); + value[i_item] = convert_result(std::move(result)); + + tokens.release(t); + }); + + return value; + } +}; + +#endif // INTERPOLATE_ITEM_VALUE_HPP -- GitLab