Skip to content
Snippets Groups Projects
Commit 1c635526 authored by Axelle Drouard's avatar Axelle Drouard
Browse files

Symmetry boundary conditions

parent bba5f847
No related branches found
No related tags found
No related merge requests found
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
#include <language/utils/InterpolateItemArray.hpp> #include <language/utils/InterpolateItemArray.hpp>
#include <mesh/Connectivity.hpp> #include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp> #include <mesh/Mesh.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp> #include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp> #include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp> #include <mesh/MeshVariant.hpp>
...@@ -228,6 +229,8 @@ class EulerKineticSolver2D ...@@ -228,6 +229,8 @@ class EulerKineticSolver2D
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr(); const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr(); const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr(); const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
const FaceValue<const TinyVector<Dimension>> nl = MeshDataManager::instance().getMeshData(m_mesh).nl();
const FaceValue<const TinyVector<Dimension>> Nl = MeshDataManager::instance().getMeshData(m_mesh).Nl();
CellValue<const double> m_inv_Vj; CellValue<const double> m_inv_Vj;
const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj(); const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
...@@ -245,7 +248,8 @@ class EulerKineticSolver2D ...@@ -245,7 +248,8 @@ class EulerKineticSolver2D
switch (bc_descriptor->type()) { switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: { case IBoundaryConditionDescriptor::Type::symmetry: {
bc_list.emplace_back( bc_list.emplace_back(
SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()))); SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break; break;
} }
case IBoundaryConditionDescriptor::Type::inflow_list: { case IBoundaryConditionDescriptor::Type::inflow_list: {
...@@ -469,9 +473,6 @@ class EulerKineticSolver2D ...@@ -469,9 +473,6 @@ class EulerKineticSolver2D
if (li_nlr > 0) { if (li_nlr > 0) {
Fr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr; Fr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
sum += li_nlr; sum += li_nlr;
// std::cout << "i_wave = " << i_wave << ", node_id = " << node_id << ", face_id = " << face_id << ",
// cell_id = " << cell_id << ", lambda = " << m_lambda_vector[i_wave] << ", nlr = " << n_face << ", njr
// = " << n_node << "\n";
} }
} }
if (sum != 0) { if (sum != 0) {
...@@ -484,6 +485,139 @@ class EulerKineticSolver2D ...@@ -484,6 +485,139 @@ class EulerKineticSolver2D
return Fr; return Fr;
} }
FaceArray<double>
compute_Flux1_face(const DiscreteFunctionP0Vector<const double> Fn)
{
const size_t nb_waves = m_lambda_vector.size();
FaceArray<double> Fl(m_mesh.connectivity(), nb_waves);
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
parallel_for(
p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
Fl[face_id][i_wave] = 0;
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nl[face_id];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nl > 0) {
Fl[face_id][i_wave] += Fn[cell_id][i_wave];
}
}
}
});
return Fl;
}
void
apply_scalar_bc(FaceArray<double> Fl_rho,
FaceArray<double> Fl_rho_u1,
FaceArray<double> Fl_rho_u2,
FaceArray<double> Fl_rho_E,
const DiscreteFunctionP0<const double> rho,
const DiscreteFunctionP0<const double> rho_u1,
const DiscreteFunctionP0<const double> rho_u2,
const DiscreteFunctionP0<const double> rho_E,
const BoundaryConditionList& bc_list)
{
const size_t nb_waves = m_lambda_vector.size();
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t i = 0; i < nb_waves; ++i) {
for (size_t d = 0; d < MeshType::Dimension; ++d) {
inv_S[d] += std::pow(m_lambda_vector[i][d], 2);
}
}
for (size_t d = 0; d < MeshType::Dimension; ++d) {
inv_S[d] = 1. / inv_S[d];
}
for (auto&& bc_v : bc_list) {
std::visit(
[=, this](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto face_list = bc.faceList();
auto n = bc.outgoingNormal();
auto nxn = tensorProduct(n, n);
TinyMatrix<Dimension> I = identity;
auto txt = I - nxn;
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nl[face_id];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nl < 0) {
double rhoj_prime = rho[cell_id];
TinyVector<Dimension> rho_uj{rho_u1[cell_id], rho_u2[cell_id]};
TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj;
double rho_Ej_prime = rho_E[cell_id];
double rho_e =
rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) *
(rho_u1[cell_id] * rho_u1[cell_id] + rho_u2[cell_id] * rho_u2[cell_id]);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
TinyVector<Dimension> A_rho_u1_prime{rho_uj_prime[0] * rho_uj_prime[0] / rhoj_prime + p,
rho_uj_prime[0] * rho_uj_prime[1] / rhoj_prime};
TinyVector<Dimension> A_rho_u2_prime{rho_uj_prime[0] * rho_uj_prime[1] / rhoj_prime,
rho_uj_prime[1] * rho_uj_prime[1] / rhoj_prime + p};
TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
double Fn_rho_prime = rhoj_prime / nb_waves +
inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_prime[0] +
inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_prime[1];
double Fn_rho_u1_prime = rho_uj_prime[0] / nb_waves +
inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_u1_prime[0] +
inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_u1_prime[1];
double Fn_rho_u2_prime = rho_uj_prime[1] / nb_waves +
inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_u2_prime[0] +
inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_u2_prime[1];
double Fn_rho_E_prime = rho_Ej_prime / nb_waves +
inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_E_prime[0] +
inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_E_prime[1];
Fl_rho[face_id][i_wave] += Fn_rho_prime;
Fl_rho_u1[face_id][i_wave] += Fn_rho_u1_prime;
Fl_rho_u2[face_id][i_wave] += Fn_rho_u2_prime;
Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
}
}
}
}
}
},
bc_v);
}
}
DiscreteFunctionP0Vector<double> DiscreteFunctionP0Vector<double>
compute_deltaLFn(NodeArray<double> Fr) compute_deltaLFn(NodeArray<double> Fr)
{ {
...@@ -558,6 +692,36 @@ class EulerKineticSolver2D ...@@ -558,6 +692,36 @@ class EulerKineticSolver2D
return deltaLFn; return deltaLFn;
} }
DiscreteFunctionP0Vector<double>
compute_deltaLFn_face(FaceArray<double> Fl)
{
const size_t nb_waves = m_lambda_vector.size();
DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
deltaLFn[cell_id][i_wave] = 0;
for (size_t i_face_cell = 0; i_face_cell < cell_to_face.size(); ++i_face_cell) {
const FaceId face_id = cell_to_face[i_face_cell];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
const double li_Nl = sign * dot(m_lambda_vector[i_wave], Nl[face_id]);
deltaLFn[cell_id][i_wave] += Fl[face_id][i_wave] * li_Nl;
}
}
});
return deltaLFn;
}
CellValue<bool> CellValue<bool>
getBoundaryCells(const BoundaryConditionList& bc_list) getBoundaryCells(const BoundaryConditionList& bc_list)
{ {
...@@ -668,20 +832,32 @@ class EulerKineticSolver2D ...@@ -668,20 +832,32 @@ class EulerKineticSolver2D
// NodeArray<double> Fr_rho_u2 = compute_Flux1(Fn_rho_u2); // NodeArray<double> Fr_rho_u2 = compute_Flux1(Fn_rho_u2);
// NodeArray<double> Fr_rho_E = compute_Flux1(Fn_rho_E); // NodeArray<double> Fr_rho_E = compute_Flux1(Fn_rho_E);
FaceArrayPerNode<double> Fr_rho = compute_Flux1_eucchlyd(Fn_rho); // FaceArrayPerNode<double> Fr_rho = compute_Flux1_eucchlyd(Fn_rho);
FaceArrayPerNode<double> Fr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1); // FaceArrayPerNode<double> Fr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1);
FaceArrayPerNode<double> Fr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2); // FaceArrayPerNode<double> Fr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2);
FaceArrayPerNode<double> Fr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E); // FaceArrayPerNode<double> Fr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E);
FaceArray<double> Fr_rho = compute_Flux1_face(Fn_rho);
FaceArray<double> Fr_rho_u1 = compute_Flux1_face(Fn_rho_u1);
FaceArray<double> Fr_rho_u2 = compute_Flux1_face(Fn_rho_u2);
FaceArray<double> Fr_rho_E = compute_Flux1_face(Fn_rho_E);
apply_scalar_bc(Fr_rho, Fr_rho_u1, Fr_rho_u2, Fr_rho_E, rho, rho_u1, rho_u2, rho_E, bc_list);
// DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn(Fr_rho); // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn(Fr_rho);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn(Fr_rho_u1); // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn(Fr_rho_u1);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn(Fr_rho_u2); // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn(Fr_rho_u2);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn(Fr_rho_E); // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn(Fr_rho_E);
DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Fr_rho); // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Fr_rho);
DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Fr_rho_u1); // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Fr_rho_u1);
DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Fr_rho_u2); // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Fr_rho_u2);
DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Fr_rho_E); // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Fr_rho_E);
DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_face(Fr_rho);
DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_face(Fr_rho_u1);
DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_face(Fr_rho_u2);
DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_face(Fr_rho_E);
const CellArray<const TinyVector<Dimension>> APFn = getA(rho, rho_u1, rho_u2, rho_E); const CellArray<const TinyVector<Dimension>> APFn = getA(rho, rho_u1, rho_u2, rho_E);
const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho{p_mesh}; const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho{p_mesh};
...@@ -866,6 +1042,7 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition ...@@ -866,6 +1042,7 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
private: private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public: public:
const Rd& const Rd&
...@@ -886,8 +1063,21 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition ...@@ -886,8 +1063,21 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
return m_mesh_flat_node_boundary.nodeList(); return m_mesh_flat_node_boundary.nodeList();
} }
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary) size_t
: m_mesh_flat_node_boundary(mesh_flat_node_boundary) numberOfFaces() const
{
return m_mesh_flat_face_boundary.faceList().size();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_flat_face_boundary.faceList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
{ {
; ;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment