diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 4447fdc4f22fdee97a8ae24dad75410c1aef0051..e0bd1af64b6b657686b9168dd9eb05b2469ebbfe 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -38,6 +38,7 @@
#include <scheme/LocalDtAcousticSolver.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>
#include <memory>
@@ -266,6 +267,15 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("coupling", std::function(
+
+ [](std::shared_ptr<const IBoundaryDescriptor> boundary)
+ -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+ return std::make_shared<CouplingBoundaryConditionDescriptor>(boundary);
+ }
+
+ ));
+
this->_addBuiltinFunction("pressure", std::function(
[](std::shared_ptr<const IBoundaryDescriptor> boundary,
@@ -398,21 +408,33 @@ SchemeModule::SchemeModule()
this->_addBuiltinFunction("local_dt_glace_solver",
std::function(
- [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& E,
- const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list,
- const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>> {
- return LocalDtAcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list1,
+ const double& dt1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list2,
+ const size_t& q) -> std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),getCommonMesh({rho2, u2, E2, c2, p2})}
.solver()
- .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt, rho, u, E, c, p,
- bc_descriptor_list);
+ .apply(LocalDtAcousticSolverHandler::SolverType::Glace, dt1, q, rho1, rho2, u1, u2, E1, E2, c1, c2, p1, p2,
+ bc_descriptor_list1, bc_descriptor_list2);
}
));
@@ -420,21 +442,33 @@ SchemeModule::SchemeModule()
this->_addBuiltinFunction("local_dt_eucclhyd_solver",
std::function(
- [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& E,
- const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
- bc_descriptor_list,
- const double& dt) -> std::tuple<std::shared_ptr<const IMesh>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>> {
- return LocalDtAcousticSolverHandler{getCommonMesh({rho, u, E, c, p})}
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list1,
+ const double& dt1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list2,
+ const size_t& q) -> std::tuple<std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return LocalDtAcousticSolverHandler{getCommonMesh({rho1, u1, E1, c1, p1}),getCommonMesh({rho2, u2, E2, c2, p2})}
.solver()
- .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt, rho, u, E, c, p,
- bc_descriptor_list);
+ .apply(LocalDtAcousticSolverHandler::SolverType::Eucclhyd, dt1, q, rho1, rho2, u1, u2, E1, E2, c1, c2, p1, p2,
+ bc_descriptor_list1, bc_descriptor_list2);
}
));
diff --git a/src/scheme/CouplingBoundaryConditionDescriptor.hpp b/src/scheme/CouplingBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3b4b26d320ab5318efee7ad95b00d34b53e78c67
--- /dev/null
+++ b/src/scheme/CouplingBoundaryConditionDescriptor.hpp
@@ -0,0 +1,46 @@
+#ifndef COUPLING_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define COUPLING_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#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;
+};
+
+#endif // COUPLING_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index f99eff9355f6fdf6fbd65ad1898f3cba655af284..5092f1973a9988be5342c48ebbdf0445ecb9d5b4 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -17,7 +17,8 @@ class IBoundaryConditionDescriptor
fixed,
free,
neumann,
- symmetry
+ symmetry,
+ coupling
};
protected:
diff --git a/src/scheme/LocalDtAcousticSolver.cpp b/src/scheme/LocalDtAcousticSolver.cpp
index ade62e0e25736800a9b6f4bd56c7d8a2e3f88cf2..e535c6e8dae0eddb2ebbc7d06ebe7d36cb5e4f43 100644
--- a/src/scheme/LocalDtAcousticSolver.cpp
+++ b/src/scheme/LocalDtAcousticSolver.cpp
@@ -15,6 +15,7 @@
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <utils/Socket.hpp>
#include <variant>
@@ -39,12 +40,14 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
class ExternalFSIVelocityBoundaryCondition;
+ class CouplingBoundaryCondition;
using BoundaryCondition = std::variant<FixedBoundaryCondition,
PressureBoundaryCondition,
SymmetryBoundaryCondition,
VelocityBoundaryCondition,
- ExternalFSIVelocityBoundaryCondition>;
+ ExternalFSIVelocityBoundaryCondition,
+ CouplingBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -268,6 +271,11 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
}
break;
}
+ case IBoundaryConditionDescriptor::Type::coupling: {
+ bc_list.emplace_back(
+ CouplingBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
default: {
is_valid_boundary_condition = false;
}
@@ -289,6 +297,20 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
const DiscreteScalarFunction& p,
NodeValue<Rdxd>& Ar,
NodeValue<Rd>& br) const;
+ void _applyCouplingBC(const BoundaryConditionList& bc_list1,
+ const BoundaryConditionList& bc_list2,
+ const MeshType& mesh1,
+ const MeshType& mesh2,
+ NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& br1,
+ NodeValue<Rd>& br2) const;
+ void _applyCouplingBC(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& ur,
+ NodeValue<Rd>& CR_ur,
+ NodeValuePerCell<Rd>& Fjr,
+ NodeValuePerCell<Rd>& CR_Fjr) const;
void
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
@@ -301,7 +323,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
this->_applyVelocityBC(bc_list, Ar, br);
}
- NodeValue<const Rd>
+ NodeValue<Rd>
_computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
{
NodeValue<Rd> u{mesh.connectivity()};
@@ -380,10 +402,140 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
}
+ std::tuple<const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ const std::shared_ptr<const ItemValueVariant>,
+ const std::shared_ptr<const SubItemValuePerItemVariant>,
+ NodeValue<Rd>,
+ NodeValuePerCell<Rd>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
+ {
+ std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, c1_v, u1_v, p1_v});
+ std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, c2_v, u2_v, p2_v});
+ if (not i_mesh1) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho1_v, c1_v, u1_v, p1_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ if (not i_mesh2) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho2_v, c2_v, u2_v, p2_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ const MeshType& mesh1 = dynamic_cast<const MeshType&>(*i_mesh1);
+ const DiscreteScalarFunction& rho1 = rho1_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& c1 = c1_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u1 = u1_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& p1 = p1_v->get<DiscreteScalarFunction>();
+
+ const MeshType& mesh2 = dynamic_cast<const MeshType&>(*i_mesh2);
+ const DiscreteScalarFunction& rho2 = rho2_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& c2 = c2_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u2 = u2_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& p2 = p2_v->get<DiscreteScalarFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * c1);
+ NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * c2);
+
+ NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+ NodeValue<Rd> br1 = this->_computeBr(mesh1, Ajr1, u1, p1);
+ NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+ NodeValue<Rd> br2 = this->_computeBr(mesh2, Ajr2, u2, p2);
+
+ const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+ const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+ this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+ this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+ this->_applyExternalVelocityBC(bc_list1, p1, Ar1, br1);
+ this->_applyExternalVelocityBC(bc_list2, p2, Ar2, br2);
+ this->_applyCouplingBC(bc_list1,bc_list2,mesh1,mesh2,Ar1,Ar2,br1,br2);
+
+ synchronize(Ar1);
+ synchronize(br1);
+ synchronize(Ar2);
+ synchronize(br2);
+
+ NodeValue<Rd> ur1 = this->_computeUr(mesh1, Ar1, br1);
+ NodeValuePerCell<Rd> Fjr1 = this->_computeFjr(mesh1, Ajr1, ur1, u1, p1);
+ NodeValue<Rd> ur2 = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> Fjr2 = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
+ NodeValue<Rd> CR_ur = this->_computeUr(mesh2, Ar2, br2);
+ NodeValuePerCell<Rd> CR_Fjr = this->_computeFjr(mesh2, Ajr2, ur2, u2, p2);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+ std::make_shared<const ItemValueVariant>(ur2),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
+ CR_ur,
+ CR_Fjr);
+ }
+
+ std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+ compute_fluxes(const SolverType& solver_type,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ NodeValue<Rd> CR_ur,
+ NodeValuePerCell<Rd> CR_Fjr) const
+ {
+ std::shared_ptr i_mesh = getCommonMesh({rho_v, c_v, u_v, p_v});
+ if (not i_mesh) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, c_v, u_v, p_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*i_mesh);
+ const DiscreteScalarFunction& rho = rho_v->get<DiscreteScalarFunction>();
+ const DiscreteScalarFunction& c = c_v->get<DiscreteScalarFunction>();
+ const DiscreteVectorFunction& u = u_v->get<DiscreteVectorFunction>();
+ const DiscreteScalarFunction& p = p_v->get<DiscreteScalarFunction>();
+
+ NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * c);
+
+ NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+ NodeValue<Rd> br = this->_computeBr(mesh, Ajr, u, p);
+
+ const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+ this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+ this->_applyExternalVelocityBC(bc_list, p, Ar, br);
+
+ synchronize(Ar);
+ synchronize(br);
+
+ NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
+
+ this->_applyCouplingBC(bc_list,mesh,ur,CR_ur,Fjr,CR_Fjr);
+
+ return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+ std::make_shared<const SubItemValuePerItemVariant>(Fjr));
+ }
+
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const DiscreteFunctionVariant>>
apply_fluxes(const double& dt,
const MeshType& mesh,
const DiscreteScalarFunction& rho,
@@ -457,6 +609,7 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
throw NormalError("acoustic solver expects P0 functions");
}
+
return this->apply_fluxes(dt, //
dynamic_cast<const MeshType&>(*i_mesh), //
rho_v->get<DiscreteScalarFunction>(), //
@@ -469,18 +622,55 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver final
std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
- const double& dt,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& E,
- const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ const double& dt1,
+ const size_t& q,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const
{
- auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
- return apply_fluxes(dt, rho, u, E, ur, Fjr);
+ std::shared_ptr<const IMesh> new_m2 = getCommonMesh({rho2, c2, u2, p2});
+ std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_u2 = u2;
+ std::shared_ptr<const DiscreteFunctionVariant> new_E2 = E2;
+
+ auto [ur1, Fjr1, ur2, Fjr2,CR_ur,CR_Fjr] = compute_fluxes(solver_type, rho1, c1, u1, p1, bc_descriptor_list1, rho2, c2, u2, p2, bc_descriptor_list2);
+ auto [new_m1,new_rho1,new_u1,new_E1] = apply_fluxes(dt1, rho1, u1, E1, ur1, Fjr1);
+
+ const double& dt2 = dt1/q;
+ const double& gamma = 1.4;
+
+ for(size_t i = 0; i<q; i++){
+ std::tie(new_m2,new_rho2,new_u2,new_E2) = apply_fluxes(dt2,new_rho2,new_u2,new_E2,ur2,Fjr2);
+
+ DiscreteScalarFunction rho_d = new_rho2->get<DiscreteScalarFunction>();
+ DiscreteVectorFunction u_d = new_u2->get<DiscreteVectorFunction>();
+ DiscreteScalarFunction E_d = new_E2->get<DiscreteScalarFunction>();
+ DiscreteScalarFunction eps = E_d - 0.5 * dot(u_d,u_d);
+ DiscreteScalarFunction p_d = (gamma-1) * rho_d * eps;
+ DiscreteScalarFunction c_d = sqrt(gamma * p_d / rho_d);
+
+ std::shared_ptr<const DiscreteFunctionVariant> new_p2 = std::make_shared<DiscreteFunctionVariant>(p_d);
+ std::shared_ptr<const DiscreteFunctionVariant> new_c2 = std::make_shared<DiscreteFunctionVariant>(c_d);
+
+ auto [ur2,Fjr2] = compute_fluxes(solver_type,new_rho2,new_c2,new_u2,new_p2,bc_descriptor_list2,CR_ur,CR_Fjr);
+ }
+
+ return {new_m1,new_rho1,new_u1,new_E1,new_m2,new_rho2,new_u2,new_E2};
}
LocalDtAcousticSolver() = default;
@@ -651,6 +841,125 @@ LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyExternalVe
}
}
+template <size_t Dimension>
+void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list1,
+ const BoundaryConditionList& bc_list2,
+ const MeshType& mesh1,
+ const MeshType& mesh2,
+ NodeValue<Rdxd>& Ar1,
+ NodeValue<Rdxd>& Ar2,
+ NodeValue<Rd>& br1,
+ NodeValue<Rd>& br2) const
+{
+ NodeValue<Rd> xr1 = copy(mesh1.xr());
+ NodeValue<Rd> xr2 = copy(mesh2.xr());
+ for(const auto& boundary_condition1 : bc_list1){
+ std::visit([&](auto&& bc1) {
+ using T1 = std::decay_t<decltype(bc1)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition,T1>){
+ const auto& node_list1 = bc1.nodeList();
+ for(const auto& boundary_condition2 : bc_list2){
+ std::visit([&](auto&& bc2) {
+ using T2 = std::decay_t<decltype(bc2)>;
+ if constexpr (std::is_same_v<CouplingBoundaryCondition,T2>){
+ const auto& node_list2 = bc2.nodeList();
+ if constexpr (Dimension == 1){
+ for(size_t i = 0; i<node_list1.size(); i++){
+ NodeId node_id1 = node_list1[i];
+ NodeId node_id2;
+ for(size_t j = 0; j<node_list2.size(); j++){
+ node_id2 = node_list2[j];
+ if(abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-14){
+ Ar1[node_id1] += Ar2[node_id2];
+ Ar2[node_id2] = Ar1[node_id1];
+
+ br1[node_id1] += br2[node_id2];
+ br2[node_id2] = br1[node_id1];
+
+ break;
+ }
+ }
+ }
+ } else if constexpr (Dimension == 2){
+ for(size_t i = 0; i<node_list1.size(); i++){
+ NodeId node_id1 = node_list1[i];
+ NodeId node_id2;
+ for(size_t j = 0; j<node_list2.size(); j++){
+ node_id2 = node_list2[j];
+ if((abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-14) and
+ (abs(xr1[node_id1][1] - xr2[node_id2][1]) < 1e-14)){
+ Ar1[node_id1] += Ar2[node_id2];
+ Ar2[node_id2] = Ar1[node_id1];
+
+ br1[node_id1] += br2[node_id2];
+ br2[node_id2] = br1[node_id1];
+
+ break;
+ }
+ }
+ }
+ } else {
+ for(size_t i = 0; i<node_list1.size(); i++){
+ NodeId node_id1 = node_list1[i];
+ NodeId node_id2;
+ for(size_t j = 0; j<node_list2.size(); j++){
+ node_id2 = node_list2[j];
+ if((abs(xr1[node_id1][0] - xr2[node_id2][0]) < 1e-14) and
+ (abs(xr1[node_id1][1] - xr2[node_id2][1]) < 1e-14) and
+ (abs(xr1[node_id1][2] - xr2[node_id2][2]) < 1e-14)){
+ Ar1[node_id1] += Ar2[node_id2];
+ Ar2[node_id2] = Ar1[node_id1];
+
+ br1[node_id1] += br2[node_id2];
+ br2[node_id2] = br1[node_id1];
+
+ break;
+ }
+ }
+ }
+ }
+ }
+ },boundary_condition2);
+ }
+ }
+ }, boundary_condition1);
+ }
+}
+
+template <size_t Dimension>
+void LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::_applyCouplingBC(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValue<Rd>& ur,
+ NodeValue<Rd>& CR_ur,
+ NodeValuePerCell<Rd>& Fjr,
+ NodeValuePerCell<Rd>& CR_Fjr) 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();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ for(size_t i = 0; i < node_list.size(); i++){
+ NodeId node_id = node_list[i];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ for(size_t j = 0; j < node_to_cell.size(); j++){
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+ Fjr(J,R) = CR_Fjr(J,R);
+ }
+ ur[node_id] = CR_ur[node_id];
+ }
+ }
+ }, boundary_condition);
+ }
+}
+
template <size_t Dimension>
class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::FixedBoundaryCondition
{
@@ -841,13 +1150,45 @@ class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<Dimension>::SymmetryBo
~SymmetryBoundaryCondition() = default;
};
-LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh)
+template <size_t Dimension>
+class LocalDtAcousticSolverHandler::LocalDtAcousticSolver<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;
+
+};
+
+LocalDtAcousticSolverHandler::LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh1, const std::shared_ptr<const IMesh>& i_mesh2)
{
- if (not i_mesh) {
- throw NormalError("discrete functions are not defined on the same mesh");
+ if (not i_mesh1) {
+ throw NormalError("mesh1 not defined");
}
- switch (i_mesh->dimension()) {
+ if (not i_mesh2) {
+ throw NormalError("mesh2 not defined");
+ }
+
+ if (not (i_mesh1->dimension()==i_mesh2->dimension())) {
+ throw NormalError("mesh1 and mesh2 must have the same dimension");
+ }
+
+
+ switch (i_mesh1->dimension()) {
case 1: {
m_acoustic_solver = std::make_unique<LocalDtAcousticSolver<1>>();
break;
diff --git a/src/scheme/LocalDtAcousticSolver.hpp b/src/scheme/LocalDtAcousticSolver.hpp
index 101054f7974aa360c16e655b0f3b9f95781fa426..3f84ae440822e22c8be14ddcb6dc60d1858e3723 100644
--- a/src/scheme/LocalDtAcousticSolver.hpp
+++ b/src/scheme/LocalDtAcousticSolver.hpp
@@ -47,15 +47,26 @@ class LocalDtAcousticSolverHandler
virtual std::tuple<std::shared_ptr<const IMesh>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const IMesh>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
apply(const SolverType& solver_type,
- const double& dt,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho,
- const std::shared_ptr<const DiscreteFunctionVariant>& u,
- const std::shared_ptr<const DiscreteFunctionVariant>& E,
- const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+ const double& dt1,
+ const size_t& q,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p1,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p2,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2) const = 0;
ILocalDtAcousticSolver() = default;
ILocalDtAcousticSolver(ILocalDtAcousticSolver&&) = default;
@@ -76,7 +87,7 @@ class LocalDtAcousticSolverHandler
return *m_acoustic_solver;
}
- LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& mesh);
+ LocalDtAcousticSolverHandler(const std::shared_ptr<const IMesh>& mesh1,const std::shared_ptr<const IMesh>& mesh2);
};
#endif // LOCAL_DT_ACOUSTIC_SOLVER_HPP