From 90701aa9381d5ed345c1524970e0b306b1b551cb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Sun, 21 May 2023 22:23:52 +0200
Subject: [PATCH] Sets inflow boundary conditions for P0 vectors
---
src/scheme/FluxingAdvectionSolver.cpp | 101 +++++++++++++++++++++-----
1 file changed, 84 insertions(+), 17 deletions(-)
diff --git a/src/scheme/FluxingAdvectionSolver.cpp b/src/scheme/FluxingAdvectionSolver.cpp
index 97aad4bb8..264ecef04 100644
--- a/src/scheme/FluxingAdvectionSolver.cpp
+++ b/src/scheme/FluxingAdvectionSolver.cpp
@@ -1,6 +1,7 @@
#include <scheme/FluxingAdvectionSolver.hpp>
#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
#include <language/utils/InterpolateItemValue.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/IMesh.hpp>
@@ -125,9 +126,9 @@ class FluxingAdvectionSolver
break;
}
default: {
- std::ostringstream os;
- os << "invalid boundary condition for advection: " << *i_bc;
- throw NotImplementedError("inflow boundary condition");
+ std::ostringstream error_msg;
+ error_msg << "invalid boundary condition for advection: " << *i_bc;
+ throw NormalError(error_msg.str());
}
}
}
@@ -152,13 +153,26 @@ class FluxingAdvectionSolver
break;
}
case IBoundaryConditionDescriptor::Type::inflow: {
- throw NotImplementedError("inflow boundary condition");
+ const InflowBoundaryConditionDescriptor& inflow_bc_descriptor =
+ dynamic_cast<const InflowBoundaryConditionDescriptor&>(*i_bc);
+
+ MeshFaceBoundary<Dimension> mesh_face_boundary =
+ getMeshFaceBoundary(mesh, inflow_bc_descriptor.boundaryDescriptor());
+
+ FaceValue<const Rd> xl = MeshDataManager::instance().getMeshData(*m_old_mesh).xl();
+ Table<const double> array_list =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>({inflow_bc_descriptor
+ .functionSymbolId()},
+ xl, mesh_face_boundary.faceList());
+
+ m_boundary_condition_list[i_quantity].emplace_back(
+ InflowArrayBoundaryCondition{mesh_face_boundary, array_list});
break;
}
default: {
- std::ostringstream os;
- os << "invalid boundary condition for advection: " << *i_bc;
- throw NotImplementedError("inflow boundary condition");
+ std::ostringstream error_msg;
+ error_msg << "invalid boundary condition for advection: " << *i_bc;
+ throw NormalError(error_msg.str());
}
}
}
@@ -556,13 +570,9 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
const CellId cell_id = m_donnor_cell[face_id];
const double fluxing_volume = m_cycle_fluxing_volume[face_id];
if (cell_id == std::numeric_limits<CellId::base_type>::max()) {
- if constexpr (is_item_value_v<CellDataType>) {
- Assert(face_to_cell_matrix[face_id].size() == 1, "boundary face must be connected to a single cell");
- const CellId face_cell_id = face_to_cell_matrix[face_id][0];
- new_q[face_cell_id] += fluxing_volume * value_list[i_face];
- } else {
- throw UnexpectedError("invalid boundary condition for fluxing advection");
- }
+ Assert(face_to_cell_matrix[face_id].size() == 1, "boundary face must be connected to a single cell");
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ new_q[face_cell_id] += fluxing_volume * value_list[i_face];
} else {
const CellId face_cell_id = face_to_cell_matrix[face_id][0];
if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
@@ -580,7 +590,33 @@ FluxingAdvectionSolver<Dimension>::_remapOne(const CellValue<const double>& step
}
} else if constexpr (std::is_same_v<TypeOfBC, InflowArrayBoundaryCondition>) {
if constexpr (is_item_array_v<CellDataType>) {
- throw NotImplementedError("advect array");
+ auto face_list = bc.faceList();
+ auto array_list = bc.arrayList();
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+ const CellId cell_id = m_donnor_cell[face_id];
+ const double fluxing_volume = m_cycle_fluxing_volume[face_id];
+ if (cell_id == std::numeric_limits<CellId::base_type>::max()) {
+ Assert(face_to_cell_matrix[face_id].size() == 1, "boundary face must be connected to a single cell");
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ auto new_array = new_q[face_cell_id];
+ const auto& bc_array = array_list[i_face];
+
+ for (size_t i = 0; i < new_array.size(); ++i) {
+ new_array[i] += fluxing_volume * bc_array[i];
+ }
+ } else {
+ const CellId face_cell_id = face_to_cell_matrix[face_id][0];
+ if (fluxing_volume > 1E-12 * step_Vj[face_cell_id]) {
+ std::ostringstream error_msg;
+ error_msg << "invalid inflow for face " << face_id
+ << " (number=" << m_old_mesh->connectivity().faceNumber()[face_id] << ")"
+ << "face has outflow fluxing volume " << fluxing_volume
+ << " (cell volume = " << step_Vj[face_cell_id] << ")\n";
+ throw NormalError(error_msg.str());
+ }
+ }
+ }
} else {
throw UnexpectedError("invalid boundary condition for fluxing advection");
}
@@ -673,8 +709,14 @@ FluxingAdvectionSolver<Dimension>::remap(
using DiscreteFunctionT = std::decay_t<decltype(variable)>;
if constexpr (std::is_same_v<MeshType, typename DiscreteFunctionT::MeshType>) {
this->_storeValues(variable);
- using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
- this->_storeValueBCList<DataType>(i, bc_list);
+ if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
+ using DataType = std::decay_t<typename DiscreteFunctionT::data_type>;
+ this->_storeValueBCList<DataType>(i, bc_list);
+ } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
+ this->_storeArrayBCList(i, bc_list);
+ } else {
+ throw UnexpectedError("invalid discrete function type");
+ }
} else {
throw UnexpectedError("incompatible mesh types");
}
@@ -746,6 +788,31 @@ class FluxingAdvectionSolver<Dimension>::InflowValueBoundaryCondition
template <size_t Dimension>
class FluxingAdvectionSolver<Dimension>::InflowArrayBoundaryCondition
{
+ private:
+ const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
+
+ Table<const double> m_array_list;
+
+ public:
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Table<const double>&
+ arrayList() const
+ {
+ return m_array_list;
+ }
+
+ InflowArrayBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary, Table<const double> array_list)
+ : m_mesh_face_boundary(mesh_face_boundary), m_array_list(array_list)
+ {
+ ;
+ }
+
+ ~InflowArrayBoundaryCondition() = default;
};
template <size_t Dimension>
--
GitLab