From bd3cba68e663d8cd08a781f99e95cd38e3a7cfac Mon Sep 17 00:00:00 2001
From: chantrait <teddy.chantrait@cea.fr>
Date: Fri, 30 Jun 2023 12:55:45 +0200
Subject: [PATCH] first try to couple two instances of pugs
---
src/language/modules/SchemeModule.cpp | 24 ++--
src/scheme/AcousticSolver.cpp | 115 ++++++++++++++++++
src/scheme/AcousticSolver.hpp | 10 +-
.../CouplingBoundaryConditionDescriptor.hpp | 42 +++++++
src/scheme/IBoundaryConditionDescriptor.hpp | 6 +-
src/utils/Serializer.cpp | 13 +-
6 files changed, 188 insertions(+), 22 deletions(-)
create mode 100644 src/scheme/CouplingBoundaryConditionDescriptor.hpp
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 788c4ec63..f2fc4ec5d 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -9,6 +9,7 @@
#include <language/utils/BinaryOperatorProcessorBuilder.hpp>
#include <language/utils/BuiltinFunctionEmbedder.hpp>
#include <language/utils/TypeDescriptor.hpp>
+#include <memory>
#include <mesh/Connectivity.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
#include <mesh/IZoneDescriptor.hpp>
@@ -20,6 +21,7 @@
#include <mesh/NumberedBoundaryDescriptor.hpp>
#include <scheme/AcousticSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
@@ -47,8 +49,6 @@
#include <scheme/VectorDiamondScheme.hpp>
#include <utils/Socket.hpp>
-#include <memory>
-
SchemeModule::SchemeModule()
{
this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const DiscreteFunctionVariant>>);
@@ -347,7 +347,14 @@ SchemeModule::SchemeModule()
}
));
+ this->_addBuiltinFunction("coupling", std::function(
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary)
+ -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<CouplingBoundaryConditionDescriptor>(boundary);
+ }
+
+ ));
#ifdef PUGS_HAS_COSTO
this->_addBuiltinFunction("cpl_pressure",
std::function(
@@ -360,32 +367,29 @@ SchemeModule::SchemeModule()
}
));
-#endif // PUGS_HAS_COSTO
-
- this->_addBuiltinFunction("normalstress",
+ this->_addBuiltinFunction("cpl_forces",
std::function(
[](std::shared_ptr<const IBoundaryDescriptor> boundary,
const FunctionSymbolId& normal_stress_id)
-> std::shared_ptr<const IBoundaryConditionDescriptor> {
- return std::make_shared<DirichletBoundaryConditionDescriptor>("normal-stress", boundary,
+ return std::make_shared<DirichletBoundaryConditionDescriptor>("cpl_forces", boundary,
normal_stress_id);
}
));
-#ifdef PUGS_HAS_COSTO
- this->_addBuiltinFunction("cpl_forces",
+#endif // PUGS_HAS_COSTO
+ this->_addBuiltinFunction("normalstress",
std::function(
[](std::shared_ptr<const IBoundaryDescriptor> boundary,
const FunctionSymbolId& normal_stress_id)
-> std::shared_ptr<const IBoundaryConditionDescriptor> {
- return std::make_shared<DirichletBoundaryConditionDescriptor>("cpl_forces", boundary,
+ return std::make_shared<DirichletBoundaryConditionDescriptor>("normal-stress", boundary,
normal_stress_id);
}
));
-#endif // PUGS_HAS_COSTO
this->_addBuiltinFunction("velocity", std::function(
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index f61f27cb4..29f18230a 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -7,6 +7,7 @@
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
@@ -15,6 +16,7 @@
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Messenger.hpp>
#include <utils/Socket.hpp>
#include <variant>
@@ -75,11 +77,13 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
class ExternalFSIVelocityBoundaryCondition;
+ class CouplingBoundaryCondition;
using BoundaryCondition = std::variant<FixedBoundaryCondition,
PressureBoundaryCondition,
SymmetryBoundaryCondition,
VelocityBoundaryCondition,
+ CouplingBoundaryCondition,
ExternalFSIVelocityBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -304,6 +308,12 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
}
break;
}
+ case IBoundaryConditionDescriptor::Type::coupling: {
+ const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
+ dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
+ bc_list.emplace_back(CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
default: {
is_valid_boundary_condition = false;
}
@@ -321,6 +331,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
+ void _applyCouplingBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void _applyExternalVelocityBC(const BoundaryConditionList& bc_list,
const DiscreteScalarFunction& p,
NodeValue<Rdxd>& Ar,
@@ -405,6 +416,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
this->_applyExternalVelocityBC(bc_list, p, Ar, br);
+ this->_applyCouplingBC(bc_list, Ar, br);
synchronize(Ar);
synchronize(br);
@@ -686,6 +698,87 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const
}
}
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list,
+ NodeValue<Rdxd>& Ar,
+ NodeValue<Rd>& br) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+
+ /* std::cout << "\033[01;31m" */
+ /* << "un applyCoupling" */
+ /* << "Dimension: " << Dimension << "\033[00;00m" << std::endl; */
+ /* std::cout << "\033[01;31m" */
+ /* << "node_list.size()" << node_list.size() << "\033[00;00m" << std::endl; */
+
+#ifdef PUGS_HAS_COSTO
+ auto costo = parallel::Messenger::getInstance().myCoupling;
+ const int dest = costo->myGlobalSize() - 1;
+ int tag = 200;
+ std::vector<int> shape;
+ shape.resize(3);
+ shape[0] = node_list.size();
+ shape[1] = Dimension + Dimension * Dimension;
+ /* shape[1] = Dimension; */
+ shape[2] = 0;
+
+ std::vector<double> data;
+ data.resize(shape[0] * shape[1] + shape[2]);
+ Array<TinyVector<Dimension>> br_extracted(node_list.size());
+ Array<TinyMatrix<Dimension>> Ar_extracted(node_list.size());
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
+ br_extracted[i_node] = br[node_id];
+ Ar_extracted[i_node] = Ar[node_id];
+ }
+ });
+ /*TODO: serializer directement Ar et br pour eviter une copie*/
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+ data[i_node * shape[1] + i_dim] = br_extracted[i_node][i_dim];
+ for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
+ data[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim] = Ar_extracted[i_node](i_dim, j_dim);
+ }
+ }
+ }
+ costo->sendData(shape, data, dest, tag);
+ std::vector<double> recv;
+
+ tag = 210;
+ costo->recvData(recv, dest, tag);
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+ br_extracted[i_node][i_dim] = recv[i_node * shape[1] + i_dim];
+ for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
+ Ar_extracted[i_node](i_dim, j_dim) = recv[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim];
+ }
+ }
+ }
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
+ br[node_id] = br_extracted[i_node];
+ Ar[node_id] = Ar_extracted[i_node];
+ }
+ });
+
+#endif PUGS_HAS_COSTO
+ }
+ },
+ boundary_condition);
+ }
+}
+
template <size_t Dimension>
class AcousticSolverHandler::AcousticSolver<Dimension>::FixedBoundaryCondition
{
@@ -876,6 +969,28 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryConditio
~SymmetryBoundaryCondition() = default;
};
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::CouplingBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary)
+ : m_mesh_node_boundary{mesh_node_boundary}
+ {
+ ;
+ }
+
+ ~CouplingBoundaryCondition() = default;
+};
+
AcousticSolverHandler::AcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh)
{
if (not i_mesh) {
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 36b62a32a..7d5c5e590 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,6 +1,4 @@
-#ifndef ACOUSTIC_SOLVER_HPP
-#define ACOUSTIC_SOLVER_HPP
-
+#pragma once
#include <memory>
#include <tuple>
#include <vector>
@@ -59,8 +57,8 @@ class AcousticSolverHandler
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
- IAcousticSolver() = default;
- IAcousticSolver(IAcousticSolver&&) = default;
+ IAcousticSolver() = default;
+ IAcousticSolver(IAcousticSolver&&) = default;
IAcousticSolver& operator=(IAcousticSolver&&) = default;
virtual ~IAcousticSolver() = default;
@@ -80,5 +78,3 @@ class AcousticSolverHandler
AcousticSolverHandler(const std::shared_ptr<const IMesh>& mesh);
};
-
-#endif // ACOUSTIC_SOLVER_HPP
diff --git a/src/scheme/CouplingBoundaryConditionDescriptor.hpp b/src/scheme/CouplingBoundaryConditionDescriptor.hpp
new file mode 100644
index 000000000..8880dde8e
--- /dev/null
+++ b/src/scheme/CouplingBoundaryConditionDescriptor.hpp
@@ -0,0 +1,42 @@
+#pragma once
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+
+class CouplingBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+ std::ostream&
+ _write(std::ostream& os) const final
+ {
+ os << "coupling(" << *m_boundary_descriptor << ")";
+ return os;
+ }
+
+ std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+ const IBoundaryDescriptor&
+ boundaryDescriptor() const final
+ {
+ return *m_boundary_descriptor;
+ }
+
+ Type
+ type() const final
+ {
+ return Type::coupling;
+ }
+
+ CouplingBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+ : m_boundary_descriptor(boundary_descriptor)
+ {
+ ;
+ }
+
+ CouplingBoundaryConditionDescriptor(const CouplingBoundaryConditionDescriptor&) = delete;
+ CouplingBoundaryConditionDescriptor(CouplingBoundaryConditionDescriptor&&) = delete;
+
+ ~CouplingBoundaryConditionDescriptor() = default;
+};
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 1350aa4db..431f6b47f 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -1,5 +1,4 @@
-#ifndef I_BOUNDARY_CONDITION_DESCRIPTOR_HPP
-#define I_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#pragma once
#include <iostream>
@@ -11,6 +10,7 @@ class IBoundaryConditionDescriptor
enum class Type
{
axis,
+ coupling,
dirichlet,
external,
fourier,
@@ -42,5 +42,3 @@ class IBoundaryConditionDescriptor
virtual ~IBoundaryConditionDescriptor() = default;
};
-
-#endif // I_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/utils/Serializer.cpp b/src/utils/Serializer.cpp
index 4e20ab62e..dec86016c 100644
--- a/src/utils/Serializer.cpp
+++ b/src/utils/Serializer.cpp
@@ -30,6 +30,8 @@ Serializer::apply(std::shared_ptr<const IMesh> i_mesh, const int dest, const int
positions[r][i] = xr[r][i];
}
});
+
+ /*TODO: serializer directement position pour eviter une copie*/
pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
@@ -56,6 +58,8 @@ Serializer::apply(std::shared_ptr<const IMesh> i_mesh, const int dest, const int
positions[r][i] = xr[r][i];
}
});
+
+ /*TODO: serializer directement position pour eviter une copie*/
pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
@@ -83,7 +87,7 @@ Serializer::apply(std::shared_ptr<const IMesh> i_mesh, const int dest, const int
}
});
- std::vector<double> pts;
+ /*TODO: serializer directement position pour eviter une copie*/
pts.resize(mesh->numberOfNodes() * MeshType::Dimension);
for (unsigned short r = 0; r < mesh->numberOfNodes(); ++r) {
for (unsigned short j = 0; j < MeshType::Dimension; ++j) {
@@ -144,6 +148,7 @@ Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
/* } */
});
+ /*TODO: serializer directement fposition pour eviter une copie*/
MeshFaceBoundary<1> mesh_face_boundary = getMeshFaceBoundary(*mesh, *boundary);
/* mesh_face_boundary.faceList() */
const auto face_list = mesh_face_boundary.faceList();
@@ -183,6 +188,7 @@ Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
/* } */
});
+ /*TODO: serializer directement fposition pour eviter une copie*/
MeshFaceBoundary<2> mesh_face_boundary = getMeshFaceBoundary(*mesh, *boundary);
/* mesh_face_boundary.faceList() */
const auto face_list = mesh_face_boundary.faceList();
@@ -222,6 +228,7 @@ Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
/* } */
});
+ /*TODO: serializer directement fposition pour eviter une copie*/
MeshFaceBoundary<3> mesh_face_boundary = getMeshFaceBoundary(*mesh, *boundary);
/* mesh_face_boundary.faceList() */
const auto face_list = mesh_face_boundary.faceList();
@@ -266,6 +273,7 @@ Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
/* } */
});
+ /*TODO: serializer directement position pour eviter une copie*/
MeshNodeBoundary<1> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
/* mesh_face_boundary.faceList() */
const auto node_list = mesh_node_boundary.nodeList();
@@ -302,6 +310,7 @@ Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
/* } */
});
+ /*TODO: serializer directement position pour eviter une copie*/
MeshNodeBoundary<2> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
/* mesh_face_boundary.faceList() */
const auto node_list = mesh_node_boundary.nodeList();
@@ -337,6 +346,7 @@ Serializer::apply(const std::shared_ptr<const IBoundaryDescriptor>& boundary,
/* } */
});
+ /*TODO: serializer directement position pour eviter une copie*/
MeshNodeBoundary<3> mesh_node_boundary = getMeshNodeBoundary(*mesh, *boundary);
/* mesh_face_boundary.faceList() */
const auto node_list = mesh_node_boundary.nodeList();
@@ -405,6 +415,7 @@ Serializer::apply(std::shared_ptr<const IMesh> i_mesh,
/* } */
});
+ /*TODO: serializer directement positions pour eviter une copie */
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
for (unsigned short i_dim = 0; i_dim < MeshType::Dimension; ++i_dim) {
--
GitLab