diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index fcebda3dedc014d935a7e73a96785762f4cb00e1..f80201c7e424728db5663045e5b663c5e306e4d4 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/src/language/modules/SchemeModule.cpp @@ -26,6 +26,7 @@ #include <scheme/DiscreteFunctionUtils.hpp> #include <scheme/DiscreteFunctionVectorIntegrator.hpp> #include <scheme/DiscreteFunctionVectorInterpoler.hpp> +#include <scheme/ExternalBoundaryConditionDescriptor.hpp> #include <scheme/FixedBoundaryConditionDescriptor.hpp> #include <scheme/FourierBoundaryConditionDescriptor.hpp> #include <scheme/FreeBoundaryConditionDescriptor.hpp> @@ -36,6 +37,7 @@ #include <scheme/NeumannBoundaryConditionDescriptor.hpp> #include <scheme/NumberedBoundaryDescriptor.hpp> #include <scheme/SymmetryBoundaryConditionDescriptor.hpp> +#include <utils/Socket.hpp> #include <memory> @@ -273,6 +275,20 @@ SchemeModule::SchemeModule() )); + this->_addBuiltinFunction("external_fsi_velocity", + std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr< + const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>, + const std::shared_ptr<const Socket>&)>>( + + [](std::shared_ptr<const IBoundaryDescriptor> boundary, + const std::shared_ptr<const Socket>& socket) + -> std::shared_ptr<const IBoundaryConditionDescriptor> { + return std::make_shared<ExternalBoundaryConditionDescriptor>("external_fsi_velocity", + boundary, socket); + } + + )); + this->_addBuiltinFunction("glace_solver", std::make_shared<BuiltinFunctionEmbedder<std::tuple< std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>, diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp index e5171badb11fa509a020025ad38e7e7bbd96281d..098b146e55efe52ffc7297b8ebb841206d9e3d43 100644 --- a/src/scheme/AcousticSolver.cpp +++ b/src/scheme/AcousticSolver.cpp @@ -8,11 +8,13 @@ #include <scheme/DirichletBoundaryConditionDescriptor.hpp> #include <scheme/DiscreteFunctionP0.hpp> #include <scheme/DiscreteFunctionUtils.hpp> +#include <scheme/ExternalBoundaryConditionDescriptor.hpp> #include <scheme/FixedBoundaryConditionDescriptor.hpp> #include <scheme/IBoundaryConditionDescriptor.hpp> #include <scheme/IDiscreteFunction.hpp> #include <scheme/IDiscreteFunctionDescriptor.hpp> #include <scheme/SymmetryBoundaryConditionDescriptor.hpp> +#include <utils/Socket.hpp> #include <variant> #include <vector> @@ -75,9 +77,13 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler class PressureBoundaryCondition; class SymmetryBoundaryCondition; class VelocityBoundaryCondition; + class ExternalFSIVelocityBoundaryCondition; - using BoundaryCondition = std:: - variant<FixedBoundaryCondition, PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>; + using BoundaryCondition = std::variant<FixedBoundaryCondition, + PressureBoundaryCondition, + SymmetryBoundaryCondition, + VelocityBoundaryCondition, + ExternalFSIVelocityBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -254,6 +260,18 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()))); break; } + case IBoundaryConditionDescriptor::Type::external: { + if constexpr (Dimension == 1) { + const ExternalBoundaryConditionDescriptor& extern_bc_descriptor = + dynamic_cast<const ExternalBoundaryConditionDescriptor&>(*bc_descriptor); + bc_list.emplace_back( + ExternalFSIVelocityBoundaryCondition(mesh, getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + extern_bc_descriptor.socket())); + } else { + throw NotImplementedError("external FSI velocity not implemented for dimension > 1"); + } + break; + } case IBoundaryConditionDescriptor::Type::dirichlet: { const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor = dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor); @@ -320,6 +338,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler 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 _applyExternalVelocityBC(const DiscreteScalarFunction& p, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br); void _applyBoundaryConditions(const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) @@ -380,6 +399,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p); this->_applyBoundaryConditions(mesh, Ar, br); + this->_applyExternalVelocityBC(p, Ar, br); synchronize(Ar); synchronize(br); @@ -610,6 +630,34 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdx } } +template <size_t Dimension> +void +AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const DiscreteScalarFunction& p, + 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<ExternalFSIVelocityBoundaryCondition, T>) { + const auto& node_list = bc.nodeList(); + const auto& value_list = bc.valueList(p); + + 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>::FixedBoundaryCondition { @@ -685,6 +733,55 @@ class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition ~PressureBoundaryCondition() = default; }; +template <size_t Dimension> +class AcousticSolverHandler::AcousticSolver<Dimension>::ExternalFSIVelocityBoundaryCondition +{ + private: + const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix; + const MeshNodeBoundary<Dimension> m_mesh_node_boundary; + const Array<TinyVector<Dimension>> m_value_list; + const std::shared_ptr<const Socket> m_socket; + + public: + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + Array<const TinyVector<Dimension>> + valueList(const DiscreteScalarFunction& p) const + { + if (parallel::size() > 1) { + throw NotImplementedError("Parallelism is not supported"); + } + if (m_value_list.size() != 1) { + throw NotImplementedError("Non connected boundaries are not supported"); + } + + const CellId cell_id = m_node_to_cell_matrix[m_mesh_node_boundary.nodeList()[0]][0]; + + write(*m_socket, p[cell_id]); + read(*m_socket, m_value_list[0]); + return m_value_list; + } + + ExternalFSIVelocityBoundaryCondition(const Mesh<Connectivity<Dimension>>& mesh, + const MeshNodeBoundary<Dimension>& mesh_node_boundary, + const std::shared_ptr<const Socket>& socket) + : m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()}, + m_mesh_node_boundary{mesh_node_boundary}, + m_value_list(mesh_node_boundary.nodeList().size()), + m_socket{socket} + { + if constexpr (Dimension > 1) { + throw NotImplementedError("external FSI velocity bc in dimension > 1"); + } + } + + ~ExternalFSIVelocityBoundaryCondition() = default; +}; + template <size_t Dimension> class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryCondition { diff --git a/src/scheme/ExternalBoundaryConditionDescriptor.hpp b/src/scheme/ExternalBoundaryConditionDescriptor.hpp new file mode 100644 index 0000000000000000000000000000000000000000..056592931b1ce1cf6801f7455d5e8ab3f785dc9f --- /dev/null +++ b/src/scheme/ExternalBoundaryConditionDescriptor.hpp @@ -0,0 +1,62 @@ +#ifndef EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP +#define EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP + +#include <mesh/IBoundaryDescriptor.hpp> +#include <scheme/IBoundaryConditionDescriptor.hpp> +#include <utils/Socket.hpp> + +#include <memory> + +class ExternalBoundaryConditionDescriptor : public IBoundaryConditionDescriptor +{ + private: + std::ostream& + _write(std::ostream& os) const final + { + os << "external(" << m_name << ',' << *m_boundary_descriptor << ")"; + return os; + } + + const std::string_view m_name; + + std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor; + const std::shared_ptr<const Socket> m_socket; + + public: + std::string_view + name() const + { + return m_name; + } + + const std::shared_ptr<const Socket>& + socket() const + { + return m_socket; + } + + const IBoundaryDescriptor& + boundaryDescriptor() const final + { + return *m_boundary_descriptor; + } + + Type + type() const final + { + return Type::external; + } + + ExternalBoundaryConditionDescriptor(const std::string_view name, + std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor, + const std::shared_ptr<const Socket>& socket) + : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_socket{socket} + {} + + ExternalBoundaryConditionDescriptor(const ExternalBoundaryConditionDescriptor&) = delete; + ExternalBoundaryConditionDescriptor(ExternalBoundaryConditionDescriptor&&) = delete; + + ~ExternalBoundaryConditionDescriptor() = default; +}; + +#endif // EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp index c0ea161db5d5f869d8ac8517b85ae2385af1d94a..f99eff9355f6fdf6fbd65ad1898f3cba655af284 100644 --- a/src/scheme/IBoundaryConditionDescriptor.hpp +++ b/src/scheme/IBoundaryConditionDescriptor.hpp @@ -12,6 +12,7 @@ class IBoundaryConditionDescriptor { axis, dirichlet, + external, fourier, fixed, free,