diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b9b789cfa8e3d7bc5328062c6085217dbe142120..e629783ba307ec56302a5719a6eabf13f68ea006 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -25,6 +25,7 @@
#include <mesh/MeshSmoother.hpp>
#include <mesh/MeshTraits.hpp>
#include <scheme/AcousticSolver.hpp>
+#include <scheme/DDFVSolver.hpp>
#include <scheme/AxisBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionDescriptorP0.hpp>
@@ -461,6 +462,7 @@ SchemeModule::SchemeModule()
));
+
this->_addBuiltinFunction("glace_solver",
std::function(
@@ -524,6 +526,64 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("DDFV_solver",
+ std::function(
+
+ [](const double& dt,
+ const double& dt_moyen,
+ const std::shared_ptr<const DiscreteFunctionVariant>& m_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& q_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& m_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& e_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& q_r,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return DDFVSolverHandler{getCommonMesh({u_j, E_j, p_j})}
+ .solver()
+ .solve(dt, dt_moyen, m_j, p_j, u_j, E_j, q_j, m_r, p_r, u_r, e_r, q_r,
+ bc_descriptor_list);
+ }
+
+ ));
+ this->_addBuiltinFunction("init_dual_sca",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& f_j_v,
+ const std::shared_ptr<const MeshVariant> primal_mesh_v,
+ const std::shared_ptr<const MeshVariant> dual_mesh_v) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return DDFVSolverHandler{primal_mesh_v}
+ .solver()
+ .init_dual_sca(f_j_v, primal_mesh_v, dual_mesh_v);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("init_dual_vec",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& f_j_v,
+ const std::shared_ptr<const MeshVariant> primal_mesh_v,
+ const std::shared_ptr<const MeshVariant> dual_mesh_v) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ return DDFVSolverHandler{primal_mesh_v}
+ .solver()
+ .init_dual_vec(f_j_v, primal_mesh_v, dual_mesh_v);
+ }
+
+ ));
+
this->_addBuiltinFunction("apply_acoustic_fluxes",
std::function(
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index a4642f38ffbc8a07e6e341ada1ce63e54a2137cf..5e7012de488985a9798b985f75e9e7c132ef0574 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -3,6 +3,7 @@
add_library(
PugsScheme
AcousticSolver.cpp
+ DDFVSolver.cpp
DiscreteFunctionIntegrator.cpp
DiscreteFunctionInterpoler.cpp
DiscreteFunctionUtils.cpp
diff --git a/src/scheme/DDFVSolver.cpp b/src/scheme/DDFVSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..53aafc4ffaa045d753cc19588cbfa960980146b1
--- /dev/null
+++ b/src/scheme/DDFVSolver.cpp
@@ -0,0 +1,880 @@
+#include <scheme/DDFVSolver.hpp>
+
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/ItemValueVariant.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <mesh/SubItemValuePerItemVariant.hpp>
+#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/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+
+#include <variant>
+#include <vector>
+
+
+#include <algebra/CRSMatrix.hpp>
+#include <algebra/CRSMatrixDescriptor.hpp>
+#include <algebra/LinearSolver.hpp>
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <algebra/Vector.hpp>
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/DualConnectivityManager.hpp>
+#include <mesh/DualMeshManager.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/PrimalToMedianDualConnectivityDataMapper.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
+#include <utils/Array.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/PugsAssert.hpp>
+#include <mesh/PrimalToDual1DConnectivityDataMapper.hpp>
+
+#include <iostream>
+#include <rang.hpp>
+
+
+double
+ddfv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+ const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+ return std::visit(
+ [&](auto&& p_mesh) -> double {
+ const auto& mesh = *p_mesh;
+
+ using MeshType = decltype(mesh);
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ 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);
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ c.meshVariant()->variant());
+}
+
+template <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver final : public DDFVSolverHandler::IDDFVSolver
+{
+ private:
+ static constexpr size_t Dimension = MeshType::Dimension;
+
+ using Rdxd = TinyMatrix<Dimension>;
+ using Rd = TinyVector<Dimension>;
+
+ using MeshDataType = MeshData<MeshType>;
+
+ using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
+ using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
+
+ class FixedBoundaryCondition;
+ class PressureBoundaryCondition;
+ class SymmetryBoundaryCondition;
+ class VelocityBoundaryCondition;
+ class ExternalFSIVelocityBoundaryCondition;
+ class DDFVBoundaryCondition;
+
+ using BoundaryCondition = std::variant<FixedBoundaryCondition,
+ PressureBoundaryCondition,
+ SymmetryBoundaryCondition,
+ VelocityBoundaryCondition,
+ ExternalFSIVelocityBoundaryCondition,
+ DDFVBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+
+
+
+
+ BoundaryConditionList
+ _getBCList(const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ BoundaryConditionList bc_list;
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::dirichlet: {
+ const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+ dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (dirichlet_bc_descriptor.name() == "velocity") {
+ MeshNodeBoundary mesh_node_boundary = getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
+
+ Array<const Rd> value_list =
+ InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
+ mesh.xr(),
+ mesh_node_boundary.nodeList());
+
+ bc_list.emplace_back(DDFVBoundaryCondition{mesh_node_boundary, value_list});
+ } 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 ddfv solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ void _applyDDFVBC(const BoundaryConditionList& bc_list,NodeValue<CellId>& computed_node_dual_cell_id, CellValue<Rd>& br) const;
+
+ public:
+
+
+
+
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ solve(const double& dt,///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ const double& dt_moyen,
+ const MeshType& mesh_p,
+ const MeshType& mesh_d,
+ const DiscreteScalarFunction& m_j,
+ const DiscreteScalarFunction& p_j,
+ const DiscreteVectorFunction& u_j,
+ const DiscreteScalarFunction& E_j,
+ const DiscreteScalarFunction& q_j,
+ const DiscreteScalarFunction& m_r,
+ const DiscreteScalarFunction& p_r,
+ const DiscreteVectorFunction& u_r,
+ const DiscreteScalarFunction& e_r,
+ const DiscreteScalarFunction& q_r,
+ const BoundaryConditionList& bc_list) const
+ {
+ TinyMatrix<Dimension> I =identity;
+ double gamma = 1.4;
+
+
+ const auto& cell_to_node_matrix_p = mesh_p.connectivity().cellToNodeMatrix();
+ const auto& node_to_cell_matrix_p = mesh_p.connectivity().nodeToCellMatrix();
+ const auto node_local_numbers_in_their_cells = mesh_p.connectivity().nodeLocalNumbersInTheirCells();
+ //const auto& cell_to_node_matrix_d = mesh_d.connectivity().cellToNodeMatrix();
+
+ MeshDataType& mesh_data_p = MeshDataManager::instance().getMeshData(mesh_p);
+ //MeshDataType& mesh_data_d = MeshDataManager::instance().getMeshData(mesh_d);
+ CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh_p).Vj();
+ CellValue<const double> Vr = MeshDataManager::instance().getMeshData(mesh_d).Vj();
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data_p.Cjr();
+
+ auto data_mapper_d =
+ DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(mesh_d.connectivity());
+ auto data_mapper_p =
+ DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(mesh_p.connectivity());
+
+ CellValue<double> new_p_j {mesh_p.connectivity()};
+ new_p_j=copy(p_j.cellValues());
+ CellValue<Rd> new_u_r{mesh_d.connectivity()};
+ new_u_r = copy(u_r.cellValues());
+
+//on passe du sommet en duale au cellule en primale
+ NodeValue<CellId> computed_node_primal_cell_id{mesh_d.connectivity()};
+ CellValue<CellId> primal_cell_id{mesh_p.connectivity()};
+ parallel_for(
+ mesh_p.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { primal_cell_id[cell_id] = cell_id; });
+ data_mapper_d->fromDualCell(primal_cell_id, computed_node_primal_cell_id);
+//////////////////////////////////////////////////////////////////////////////////////////
+//on passe du sommet en primale au cellule en duale
+ NodeValue<CellId> computed_node_dual_cell_id{mesh_p.connectivity()};
+ CellValue<CellId> dual_cell_id{mesh_d.connectivity()};
+ parallel_for(
+ mesh_d.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
+ data_mapper_p->fromDualCell(dual_cell_id, computed_node_dual_cell_id);
+/////////////////////////////////////////////////////////////////////////////////////////
+/*
+std::cout<<"pppppppppppppp"<<std::endl;
+ parallel_for(
+ mesh_p.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+ std::cout<<"nodeId="<<r<<"cell_id="<<dual_cell_id_r<<std::endl;
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+
+ size_t i_node_in_face = node_local_number_in_their_cells[r][node_cells[j]];
+
+ new_u_r[dual_cell_id_r] += dt_moyen * new_p_j[node_cells[j]] * (1.0/m_r[dual_cell_id_r]) * Cjr(node_cells[j],i_node_in_face) ;
+ }
+ });*/
+
+ for(NodeId r = 0; r<mesh_p.numberOfNodes(); r++){
+
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+ const unsigned int R = node_local_number_in_its_cells[j]; //numérotation local je comprends rien
+
+ new_u_r[dual_cell_id_r] -= dt_moyen * (-new_p_j[node_cells[j]] - q_j[node_cells[j]]) * (1.0/m_r[dual_cell_id_r]) * Cjr(node_cells[j],R) ;
+
+ }
+ };
+
+
+ _applyDDFVBC(bc_list, computed_node_dual_cell_id, new_u_r);
+
+
+ NodeValue<Rd> new_xr = copy(mesh_p.xr());
+
+
+ parallel_for(
+ mesh_p.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+ new_xr[r] += dt * new_u_r[dual_cell_id_r];
+ });
+
+
+ std::shared_ptr<const MeshType> new_mesh_p = std::make_shared<MeshType>(mesh_p.shared_connectivity(), new_xr);
+ std::shared_ptr<const MeshVariant> new_mesh_d_v = DualMeshManager::instance().getMedianDualMesh(std::make_shared<MeshVariant>(new_mesh_p));
+ std::shared_ptr<const MeshType> new_mesh_d =new_mesh_d_v->get<MeshType>();
+
+ MeshDataType& new_mesh_data_p = MeshDataManager::instance().getMeshData(*new_mesh_p);
+ MeshDataType& new_mesh_data_d = MeshDataManager::instance().getMeshData(*new_mesh_d);
+ const NodeValuePerCell<const Rd> new_Cjr = new_mesh_data_p.Cjr();
+ const NodeValuePerCell<const Rd> new_Crj = new_mesh_data_d.Cjr();
+
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh_p).Vj();
+ CellValue<const double> new_Vr = MeshDataManager::instance().getMeshData(*new_mesh_d).Vj();
+ //CellValue<double> new_Vr {mesh_d.connectivity()};
+
+ CellValue<double> new_rho_j{mesh_p.connectivity()};
+
+ parallel_for(
+ mesh_p.numberOfCells(), PUGS_LAMBDA(CellId j) {
+
+ new_rho_j[j] = m_j[j] / new_Vj[j];
+
+ });
+
+ CellValue<double> new_rho_r{mesh_d.connectivity()};
+ /*
+ for(NodeId r = 0; r<mesh_p.numberOfNodes(); r++){
+
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ new_Vr[dual_cell_id_r] = Vr[dual_cell_id_r];
+
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+ const unsigned int R = node_local_number_in_its_cells[j]; //numérotation local je comprends rien
+
+ new_Vr[dual_cell_id_r] -= dt_moyen * dot(u_j[node_cells[j]] , Cjr(node_cells[j],R)) ;
+
+ }
+ };
+ */
+ parallel_for(
+ mesh_d.numberOfCells(), PUGS_LAMBDA(CellId r) {
+ new_rho_r[r] = m_r[r] *2.0 / (Vr[r] + new_Vr[r]);
+ });
+
+
+ CellValue<double> new_p_r{mesh_d.connectivity()};
+ CellValue<double> new_epsilon_r = copy(e_r.cellValues());
+
+ for(NodeId r = 0; r<mesh_p.numberOfNodes(); r++){
+
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId cell_r = computed_node_dual_cell_id[r];
+ new_epsilon_r[cell_r] = new_epsilon_r[cell_r] * m_r[cell_r];
+
+ double mem=m_r[cell_r];
+
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+
+ CellId cell_j = node_cells[j];
+ const unsigned int R = node_local_number_in_its_cells[j]; //numérotation local je comprends rien
+ new_epsilon_r[cell_r] -= dt_moyen * (-q_r[cell_r] - 0.5 * p_r[cell_r]) * doubleDot(I,tensorProduct(u_j[cell_j],Cjr(cell_j,R)));
+
+ mem -= dt_moyen * 0.5 * (gamma-1.0) * new_rho_r[cell_r] * doubleDot(I,tensorProduct(u_j[cell_j],Cjr(cell_j,R)));
+ }
+
+ new_epsilon_r[cell_r] = new_epsilon_r[cell_r] / mem;
+
+ new_p_r[cell_r]=(gamma - 1.0) * new_epsilon_r[cell_r] * new_rho_r[cell_r];
+ };
+
+
+
+
+ CellValue<Rd> new_u_j = copy(u_j.cellValues());
+
+ parallel_for(
+ mesh_p.numberOfCells(), PUGS_LAMBDA(CellId j) {
+
+ const auto& cell_nodes = cell_to_node_matrix_p[j];
+
+ for (size_t r = 0; r < cell_nodes.size(); ++r)
+ {
+ CellId dual_cell_id_r = computed_node_dual_cell_id[cell_nodes[r]];
+
+ new_u_j[j] += dt * 0.5 * (1.0/m_j[j]) * (-new_p_r[dual_cell_id_r] - q_r[dual_cell_id_r]) * (Cjr(j,r) + new_Cjr(j,r));
+
+ }
+
+ });
+
+ CellValue<double> new_E_j = copy(E_j.cellValues());
+ CellValue<double> new_e_j{mesh_p.connectivity()} ;
+
+ parallel_for(
+ mesh_p.numberOfCells(), PUGS_LAMBDA(CellId j) {
+
+ const auto& cell_nodes = cell_to_node_matrix_p[j];
+
+ for (size_t r = 0; r < cell_nodes.size(); ++r)
+ {
+ CellId dual_cell_id_r = computed_node_dual_cell_id[cell_nodes[r]];
+
+ new_E_j[j] += dt * 0.5 * (1.0/m_j[j]) * (-new_p_r[dual_cell_id_r] - q_r[dual_cell_id_r]) * dot(new_u_r[dual_cell_id_r],(Cjr(j,r) + new_Cjr(j,r)));
+
+ }
+
+ new_e_j[j] = new_E_j[j] - 0.5 * dot(new_u_j[j],new_u_j[j]);
+ new_p_j[j] = (gamma - 1.0) * new_e_j[j] * new_rho_j[j];
+ });
+ CellValue<double> grad_u_r{mesh_d.connectivity()};
+ CellValue<double> new_q_r{mesh_d.connectivity()};
+
+ for(NodeId r = 0; r<mesh_p.numberOfNodes(); r++){
+
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+
+ grad_u_r[dual_cell_id_r] = 0;
+
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+ const unsigned int R = node_local_number_in_its_cells[j]; //numérotation local je comprends rien
+
+ grad_u_r[dual_cell_id_r] -= 0.25 * (1.0/new_Vr[dual_cell_id_r]) * dot(new_u_j[node_cells[j]] + u_j[node_cells[j]], Cjr(node_cells[j],R) + new_Cjr(node_cells[j],R));
+
+ }
+ //double c=sqrt(gamma * new_p_r[dual_cell_id_r] / new_rho_r[dual_cell_id_r]);
+ //double lambda=new_mesh_data_d.sumOverRLjr()[dual_cell_id_r];
+ //double lambda = new_Vr[dual_cell_id_r];
+
+ new_q_r[dual_cell_id_r] = - 1.0 * new_Vr[dual_cell_id_r] * new_rho_r[dual_cell_id_r] * grad_u_r[dual_cell_id_r] * sqrt(gamma * new_p_r[dual_cell_id_r] / new_rho_r[dual_cell_id_r]);
+ //new_q_r[dual_cell_id_r] = -lambda * new_rho_r[dual_cell_id_r] * (-0.0 * c + lambda * grad_u_r[dual_cell_id_r]) * grad_u_r[dual_cell_id_r];
+ };
+
+ CellValue<double> grad_u_j{mesh_p.connectivity()};
+ CellValue<double> new_q_j{mesh_p.connectivity()};
+ for(CellId j = 0; j<mesh_p.numberOfCells(); j++){
+
+ const auto& cell_nodes = cell_to_node_matrix_p[j];
+
+ grad_u_j[j] = 0.0;
+
+ for (size_t r = 0.0; r < cell_nodes.size(); ++r)
+ {
+
+ CellId dual_cell_id_r = computed_node_dual_cell_id[cell_nodes[r]];
+ grad_u_j[j] += 0.5 * (1.0/new_Vj[j]) * dot(new_u_r[dual_cell_id_r], Cjr(j,r) + new_Cjr(j,r)) ;
+
+ }
+ //double c = sqrt(gamma * new_p_j[j] / new_rho_j[j]);
+ //double lambda = new_mesh_data_p.sumOverRLjr()[j];
+ //double lambda = new_Vj[j];
+
+ new_q_j[j] = - 1.0 * new_Vj[j] * new_rho_j[j] * grad_u_j[j] * sqrt(gamma * new_p_j[j] / new_rho_j[j]);
+ //new_q_j[j] = -lambda * new_rho_j[j] * (-0.0 * c + lambda * grad_u_j[j]) * grad_u_j[j];
+
+ };
+
+
+ return {std::make_shared<MeshVariant>(new_mesh_p),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh_p, new_rho_j)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh_p, new_u_j)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh_p, new_E_j)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh_p, new_q_j)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh_d, new_rho_r)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(new_mesh_d, new_u_r)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh_d, new_epsilon_r)),
+ std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(new_mesh_d, new_q_r))};
+ }
+
+
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ solve(const double& dt,
+ const double& dt_moyen,
+ const std::shared_ptr<const DiscreteFunctionVariant>& m_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& q_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& m_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& e_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& q_r,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ std::shared_ptr mesh_p = getCommonMesh({u_j, p_j, E_j});
+ if (not mesh_p) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ std::shared_ptr mesh_d = getCommonMesh({u_r, p_r, e_r});
+ if (not mesh_d) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ const BoundaryConditionList bc_list = this->_getBCList(*mesh_p->get<MeshType>(), bc_descriptor_list);
+
+ return this->solve(dt,
+ dt_moyen,
+ *mesh_p->get<MeshType>(),
+ *mesh_d->get<MeshType>(),
+ m_j->get<DiscreteScalarFunction>(),
+ p_j->get<DiscreteScalarFunction>(),
+ u_j->get<DiscreteVectorFunction>(),
+ E_j->get<DiscreteScalarFunction>(),
+ q_j->get<DiscreteScalarFunction>(),
+ m_r->get<DiscreteScalarFunction>(),
+ p_r->get<DiscreteScalarFunction>(),
+ u_r->get<DiscreteVectorFunction>(),
+ e_r->get<DiscreteScalarFunction>(),
+ q_r->get<DiscreteScalarFunction>(),
+ bc_list);
+ }
+
+ std::shared_ptr<const DiscreteFunctionVariant>
+ init_dual_sca(const DiscreteScalarFunction& f_j,
+ const std::shared_ptr<const MeshVariant>& primal_mesh_v,
+ const std::shared_ptr<const MeshVariant>& dual_mesh_v) const {
+//std::cout<<"début test initialisation"<<std::endl;
+ const MeshType& mesh_p= *primal_mesh_v->get<MeshType>() ;
+ const MeshType& mesh_d= *dual_mesh_v->get<MeshType>() ;
+ CellValue<double> f_r{mesh_d.connectivity()} ;
+
+ auto data_mapper_p = DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(mesh_p.connectivity());
+ const auto& node_to_cell_matrix_p = mesh_p.connectivity().nodeToCellMatrix();
+
+ NodeValue<CellId> computed_node_dual_cell_id{mesh_p.connectivity()};
+ CellValue<CellId> dual_cell_id{mesh_d.connectivity()};
+ parallel_for(
+ mesh_d.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
+ data_mapper_p->fromDualCell(dual_cell_id, computed_node_dual_cell_id);
+
+ /*parallel_for(
+ mesh_p.numberOfNodes(), PUGS_LAMBDA(NodeId r) {*/
+ for(NodeId r=0 ; r<mesh_p.numberOfNodes() ; r++){
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+
+ CellId primal_cell_id_j = node_cells[j];
+
+ f_r[dual_cell_id_r] += f_j[primal_cell_id_j] ;
+ }
+ f_r[dual_cell_id_r] = (1.0 / node_cells.size()) * f_r[dual_cell_id_r]; //pas sûre de la division (1/int)
+ //std::cout<<"u("<<dual_cell_id_r<<")= "<<f_r[dual_cell_id_r]<<std::endl;
+ }
+
+ /*});*/
+
+ //std::cout<<"fin de test";
+ return std::make_shared<DiscreteFunctionVariant>(DiscreteScalarFunction(dual_mesh_v, f_r));
+ }
+
+ std::shared_ptr<const DiscreteFunctionVariant>
+ init_dual_sca(const std::shared_ptr<const DiscreteFunctionVariant>& f_j_v,
+ const std::shared_ptr<const MeshVariant>& primal_mesh_v,
+ const std::shared_ptr<const MeshVariant>& dual_mesh_v) const {
+
+ const DiscreteScalarFunction f_j = f_j_v->get<DiscreteScalarFunction>();
+
+ return this->init_dual_sca(f_j, primal_mesh_v, dual_mesh_v);
+ }
+
+ std::shared_ptr<const DiscreteFunctionVariant>
+ init_dual_vec(const DiscreteVectorFunction& f_j,
+ const std::shared_ptr<const MeshVariant>& primal_mesh_v,
+ const std::shared_ptr<const MeshVariant>& dual_mesh_v) const{
+
+ //std::cout<<"début test initialisation"<<std::endl;
+ const MeshType& mesh_p= *primal_mesh_v->get<MeshType>() ;
+ const MeshType& mesh_d= *dual_mesh_v->get<MeshType>() ;
+
+ CellValue<Rd> f_r{mesh_d.connectivity()};
+
+ auto data_mapper_p = DualConnectivityManager::instance().getPrimalToMedianDualConnectivityDataMapper(mesh_p.connectivity());
+ const auto& node_to_cell_matrix_p = mesh_p.connectivity().nodeToCellMatrix();
+
+ NodeValue<CellId> computed_node_dual_cell_id{mesh_p.connectivity()};
+ CellValue<CellId> dual_cell_id{mesh_d.connectivity()};
+ parallel_for(
+ mesh_d.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { dual_cell_id[cell_id] = cell_id; });
+ data_mapper_p->fromDualCell(dual_cell_id, computed_node_dual_cell_id);
+/*
+ parallel_for(
+ mesh_p.numberOfNodes(), PUGS_LAMBDA(NodeId r) {*/
+ for(NodeId r=0 ; r<mesh_p.numberOfNodes() ; r++){
+ const auto& node_cells = node_to_cell_matrix_p[r];
+ CellId dual_cell_id_r = computed_node_dual_cell_id[r];
+
+ for (size_t j = 0; j < node_cells.size(); ++j)
+ {
+
+ CellId primal_cell_id_j = node_cells[j];
+
+ f_r[dual_cell_id_r] += f_j[primal_cell_id_j] ;
+ }
+ f_r[dual_cell_id_r] = (1.0 / node_cells.size()) * f_r[dual_cell_id_r]; //pas sûre de la division (1/int)
+ //std::cout<<"u("<<dual_cell_id_r<<")= "<<f_r[dual_cell_id_r]<<std::endl;
+ }
+ /* });*/
+
+ //std::cout<<"fin de test";
+
+
+ return std::make_shared<DiscreteFunctionVariant>(DiscreteVectorFunction(dual_mesh_v, f_r));
+ }
+
+ std::shared_ptr<const DiscreteFunctionVariant>
+ init_dual_vec(const std::shared_ptr<const DiscreteFunctionVariant>& f_j_v,
+ const std::shared_ptr<const MeshVariant>& primal_mesh_v,
+ const std::shared_ptr<const MeshVariant>& dual_mesh_v) const {
+
+
+ const DiscreteVectorFunction f_j = f_j_v->get<DiscreteVectorFunction>();
+
+ return this->init_dual_vec(f_j, primal_mesh_v, dual_mesh_v);
+ }
+
+ DDFVSolver() = default;
+ DDFVSolver(DDFVSolver&&) = default;
+ ~DDFVSolver() = default;
+};
+
+template <MeshConcept MeshType>
+void
+DDFVSolverHandler::DDFVSolver<MeshType>::_applyDDFVBC(const BoundaryConditionList& bc_list,
+ NodeValue<CellId>& computed_node_dual_cell_id,
+ CellValue<Rd>& u_r) 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<DDFVBoundaryCondition, 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];
+ CellId node_dual_cell_id = computed_node_dual_cell_id[node_id];
+
+ u_r[node_dual_cell_id] = value;
+
+ });
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver<MeshType>::FixedBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+ ~FixedBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver<MeshType>::PressureBoundaryCondition
+{
+ private:
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Array<const double> m_value_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const MeshFaceBoundary& mesh_face_boundary, const Array<const double>& value_list)
+ : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <>
+class DDFVSolverHandler::DDFVSolver<Mesh<1>>::PressureBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const Array<const double> m_value_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const double>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ PressureBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const Array<const double>& value_list)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+ {}
+
+ ~PressureBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver<MeshType>::ExternalFSIVelocityBoundaryCondition
+{
+ private:
+ const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
+ const MeshNodeBoundary 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<Dimension>& mesh,
+ const MeshNodeBoundary& 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 <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver<MeshType>::VelocityBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+
+ const Array<const TinyVector<Dimension>> m_value_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ VelocityBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+ {}
+
+ ~VelocityBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> 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<MeshType>& mesh_flat_node_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class DDFVSolverHandler::DDFVSolver<MeshType>::DDFVBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+
+ const Array<const TinyVector<Dimension>> m_value_list;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const TinyVector<Dimension>>&
+ valueList() const
+ {
+ return m_value_list;
+ }
+
+ DDFVBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const Array<const TinyVector<Dimension>>& value_list)
+ : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
+ {}
+
+ ~DDFVBoundaryCondition() = default;
+};
+
+DDFVSolverHandler::DDFVSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v)
+{
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ std::visit(
+ [&](auto&& mesh) {
+ using MeshType = mesh_type_t<decltype(mesh)>;
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ m_ddfv_solver = std::make_unique<DDFVSolver<MeshType>>();
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/DDFVSolver.hpp b/src/scheme/DDFVSolver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8ff34bf02a1598b7e81bdfdbbef61cf25ef5e58d
--- /dev/null
+++ b/src/scheme/DDFVSolver.hpp
@@ -0,0 +1,85 @@
+#ifndef DDFV_SOLVER_HPP
+#define DDFV_SOLVER_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+class DiscreteFunctionVariant;
+class IBoundaryConditionDescriptor;
+class MeshVariant;
+class ItemValueVariant;
+class SubItemValuePerItemVariant;
+
+double ddfv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c);
+
+class DDFVSolverHandler
+{
+ public:
+ enum class SolverType
+ {
+ Glace,
+ Eucclhyd
+ };
+
+ private:
+ struct IDDFVSolver
+ {
+ virtual std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ solve(const double& dt,
+ const double& dt_moyen,
+ const std::shared_ptr<const DiscreteFunctionVariant>& m_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& q_j,
+ const std::shared_ptr<const DiscreteFunctionVariant>& m_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& e_r,
+ const std::shared_ptr<const DiscreteFunctionVariant>& q_r,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const =0;
+
+ virtual std::shared_ptr<const DiscreteFunctionVariant>
+ init_dual_sca(const std::shared_ptr<const DiscreteFunctionVariant>& f_j_v,
+ const std::shared_ptr<const MeshVariant>& primal_mesh_v,
+ const std::shared_ptr<const MeshVariant>& dual_mesh_v) const =0;
+
+ virtual std::shared_ptr<const DiscreteFunctionVariant>
+ init_dual_vec(const std::shared_ptr<const DiscreteFunctionVariant>& f_j_v,
+ const std::shared_ptr<const MeshVariant>& primal_mesh_v,
+ const std::shared_ptr<const MeshVariant>& dual_mesh_v) const =0;
+
+ IDDFVSolver() = default;
+ IDDFVSolver(IDDFVSolver&&) = default;
+ IDDFVSolver& operator=(IDDFVSolver&&) = default;
+
+ virtual ~IDDFVSolver() = default;
+ };
+
+ template <MeshConcept MeshType>
+ class DDFVSolver;
+
+ std::unique_ptr<IDDFVSolver> m_ddfv_solver;
+
+ public:
+ const IDDFVSolver&
+ solver() const
+ {
+ return *m_ddfv_solver;
+ }
+
+ DDFVSolverHandler(const std::shared_ptr<const MeshVariant>& mesh_v);
+};
+
+#endif // DDFV_SOLVER_HPP