diff --git a/src/language/algorithms/AcousticSolverAlgorithm.hpp b/src/language/algorithms/AcousticSolverAlgorithm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3fca4703f6d69c7790159519f6728479727183a5
--- /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 477e3cc80cf439ae369cb1bd5074e53116eda91b..b38373e672a736a8fe4c64d2d6125b265e1cb7a1 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 0000000000000000000000000000000000000000..16a87b4952c7cb0c2907833705e8142d9c0a7178
--- /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