27namespace HeatTransportBHE
32 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
33 std::vector<std::unique_ptr<ParameterLib::ParameterBase>>
const& parameters,
34 unsigned const integration_order,
35 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
40 :
Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
41 integration_order, std::move(process_variables),
42 std::move(secondary_variables)),
43 _process_data(std::move(process_data)),
44 _bheMeshData(std::move(bhe_mesh_data)),
45 _asm_mat_cache{_process_data._is_linear, true }
51 "The number of the given BHE properties ({:d}) are not consistent "
52 "with the number of BHE groups in the mesh ({:d}).",
53 _process_data._vec_BHE_property.size(),
54 _bheMeshData.BHE_mat_IDs.size());
58 if (material_ids ==
nullptr)
60 OGS_FATAL(
"Not able to get material IDs! ");
87 std::vector<std::vector<MeshLib::Element*>
const*> vec_var_elements;
90 std::vector<int> vec_n_components{
102 for (
int i = 0; i < n_BHEs; i++)
104 auto const number_of_unknowns =
105 visit([](
auto const& bhe) {
return bhe.number_of_unknowns; },
112 std::make_unique<MeshLib::MeshSubset const>(
_mesh, bhe_nodes));
114 std::generate_n(std::back_inserter(all_mesh_subsets),
121 vec_n_components.push_back(number_of_unknowns);
122 vec_var_elements.push_back(&bhe_elements);
126 std::make_unique<NumLib::LocalToGlobalIndexMap>(
127 std::move(all_mesh_subsets),
139 unsigned const integration_order)
142 std::unordered_map<std::size_t, BHE::BHETypes*> element_to_bhe_map;
144 for (
int i = 0; i < n_BHEs; i++)
147 for (
auto const& e : bhe_elements)
149 element_to_bhe_map[e->getID()] =
178 ranges::to<std::vector>();
188 ranges::to<std::vector>();
193 std::vector<std::size_t>
const bhes_node_ids =
195 ranges::views::transform([](
auto const*
const node)
196 {
return node->getID(); }) |
197 ranges::to<std::vector>;
206 for (
auto const id : es.getSearchedElementIDs())
214 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
215 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
218 DBUG(
"Assemble HeatTransportBHE process.");
220 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
228 "Domain Deactivation is currently not implemnted with "
229 "linear optimization.");
247 process_id, &M, &K, &b);
255 x_prev, process_id, &M, &K, &b);
268 const double t,
double const dt, std::vector<GlobalVector*>
const& x,
269 std::vector<GlobalVector*>
const& x_prev,
int const process_id,
272 DBUG(
"AssembleWithJacobian HeatTransportBHE process.");
274 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_table = {
281 process_id, &b, &Jac);
285 double const t,
double const dt, std::vector<GlobalVector*>
const& x,
288 DBUG(
"Compute heat flux for HeatTransportBHE process.");
290 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
291 dof_tables.reserve(x.size());
292 std::generate_n(std::back_inserter(dof_tables), x.size(),
293 [&]() { return _local_to_global_index_map.get(); });
310 auto const Tout_nodes_id =
312 const std::size_t n_bc_nodes = Tout_nodes_id.size();
314 for (std::size_t i = 0; i < n_bc_nodes; i++)
327 DBUG(
"Method `tespySolver' not overridden in Python script.");
331 for (std::size_t i = 0; i < n_bc_nodes; i++)
334 std::get<2>(tespy_result)[i];
336 std::get<3>(tespy_result)[i];
338 auto const tespy_has_converged = std::get<1>(tespy_result);
339 if (tespy_has_converged ==
true)
346 std::vector<GlobalVector*>
const& x,
const double t,
const double dt,
347 int const process_id)
355 auto& [time, Tin_value, Tout_value, Tout_nodes_ids, flowrate] =
363 auto const& solution = *x[process_id];
366 const std::size_t n_bc_nodes = Tout_nodes_ids.size();
367 for (std::size_t i = 0; i < n_bc_nodes; i++)
370 Tout_value[i] = solution[Tout_nodes_ids[i]];
374 auto const server_communication_result =
376 t, dt, Tin_value, Tout_value, flowrate);
380 DBUG(
"Method `serverCommunication' not overridden in Python script.");
383 auto const& [server_communication_Tin_value,
384 server_communication_flowrate] = server_communication_result;
386 std::copy(begin(server_communication_Tin_value),
387 end(server_communication_Tin_value),
389 std::copy(begin(server_communication_flowrate),
390 end(server_communication_flowrate),
395 std::vector<GlobalVector*>
const& x,
396 std::vector<GlobalVector*>
const& ,
const double t,
397 const double dt,
int const process_id)
405 auto& [time, Tin_value, Tout_value, Tout_nodes_ids, flowrate] =
413 auto const& solution = *x[process_id];
416 const std::size_t n_bc_nodes = Tout_nodes_ids.size();
417 for (std::size_t i = 0; i < n_bc_nodes; i++)
420 Tout_value[i] = solution[Tout_nodes_ids[i]];
425 t, dt, Tin_value, Tout_value, flowrate);
429 DBUG(
"Method `serverCommunication' not overridden in Python script.");
434 [[maybe_unused]]
const double t,
double const ,
435 [[maybe_unused]] std::vector<GlobalVector*>
const& x,
436 std::vector<GlobalVector*>
const& ,
int const ,
441 auto M_normal = M.getRawMatrix();
442 auto K_normal = K.getRawMatrix();
443 auto n_original_rows = K_normal.rows();
451 M_normal.conservativeResize(
452 M_normal.rows() + n_BHE_bottom_pairs + n_BHE_top_pairs,
454 K_normal.conservativeResize(
455 K_normal.rows() + n_BHE_bottom_pairs + n_BHE_top_pairs,
458 for (std::size_t i = 0; i < n_BHE_bottom_pairs; i++)
460 Eigen::SparseVector<double> M_Plus(M_normal.cols());
462 M_normal.row(n_original_rows + i) = M_Plus;
464 Eigen::SparseVector<double> K_Plus(K_normal.cols());
467 auto const [bhe_idx, first_BHE_bottom_index, second_BHE_bottom_index] =
470 K_Plus.insert(first_BHE_bottom_index) = w_val;
471 K_Plus.insert(second_BHE_bottom_index) = -w_val;
473 K_normal.row(n_original_rows + i) = K_Plus;
476 auto b_normal = b.getRawVector();
477 Eigen::SparseVector<double> b_Plus(b_normal.rows() + n_BHE_bottom_pairs +
482 for (
int i = 0; i < b_normal.innerSize(); ++i)
484 b_Plus.insert(i) = b_normal.coeff(i);
487 for (std::size_t i = 0; i < n_BHE_top_pairs; i++)
489 Eigen::SparseVector<double> M_Plus(M_normal.cols());
491 M_normal.row(n_original_rows + n_BHE_bottom_pairs + i) = M_Plus;
493 Eigen::SparseVector<double> K_Plus(K_normal.cols());
496 auto const [bhe_idx, first_BHE_top_index, second_BHE_top_index] =
499 auto first_BHE_top_index_pair = first_BHE_top_index;
500 auto second_BHE_top_index_pair = second_BHE_top_index;
502 K_Plus.insert(first_BHE_top_index_pair) =
504 K_Plus.insert(second_BHE_top_index_pair) =
507 K_normal.row(n_original_rows + n_BHE_bottom_pairs + i) = K_Plus;
510 double const T_out = (*x[0])[second_BHE_top_index_pair];
512 auto calculate_delta_T = [&](
auto& bhe)
514 auto const T_in = bhe.updateFlowRateAndTemperature(T_out, t);
517 auto delta_T = std::visit(calculate_delta_T,
520 b_Plus.insert(n_original_rows + n_BHE_bottom_pairs + i) =
524 M.getRawMatrix() = M_normal;
525 K.getRawMatrix() = K_normal;
526 b.getRawVector() = b_Plus;
529 "The Algebraic Boundary Condition is not implemented for use with "
530 "PETsc Library! Simulation will be terminated.");
535 std::vector<std::vector<MeshLib::Node*>>
const& all_bhe_nodes)
537 const int process_id = 0;
543 for (std::size_t bhe_i = 0; bhe_i < n_BHEs; bhe_i++)
545 auto const& bhe_nodes = all_bhe_nodes[bhe_i];
549 const int variable_id = bhe_i + 1;
551 std::vector<MeshLib::Node*> bhe_boundary_nodes;
555 for (
auto const& bhe_node : bhe_nodes)
558 auto const& connected_elements =
560 const std::size_t n_line_elements = std::count_if(
561 connected_elements.begin(), connected_elements.end(),
563 { return (elem->getDimension() == 1); });
565 if (n_line_elements == 1)
567 bhe_boundary_nodes.push_back(bhe_node);
571 if (bhe_boundary_nodes.size() != 2)
574 "Error!!! The BHE boundary nodes are not correctly found, "
575 "for every single BHE, there should be 2 boundary nodes.");
590 if ((*bhe_boundary_nodes[0])[2] == (*bhe_boundary_nodes[1])[2])
593 "For 1P-type BHE, the BHE inflow and outflow "
594 "nodes are identified according to their mesh node id in "
602 if ((*bhe_boundary_nodes[0])[2] < (*bhe_boundary_nodes[1])[2])
604 std::swap(bhe_boundary_nodes[0], bhe_boundary_nodes[1]);
608 auto get_global_index =
609 [&](std::size_t
const node_id,
int const component)
613 variable_id, component);
616 auto get_global_bhe_bc_indices =
618 std::pair<std::size_t ,
int >, 2>
619 nodes_and_components)
621 return std::make_pair(
622 get_global_index(nodes_and_components[0].first,
623 nodes_and_components[0].second),
624 get_global_index(nodes_and_components[1].first,
625 nodes_and_components[1].second));
628 auto get_global_bhe_bc_indices_with_bhe_idx =
629 [&](std::size_t bhe_idx,
631 std::pair<std::size_t ,
int >, 2>
632 nodes_and_components)
634 return std::make_tuple(
636 get_global_index(nodes_and_components[0].first,
637 nodes_and_components[0].second),
638 get_global_index(nodes_and_components[1].first,
639 nodes_and_components[1].second));
643 [&, bc_top_node_id = bhe_boundary_nodes[0]->getID(),
644 bc_bottom_node_id = bhe_boundary_nodes[1]->getID()](
auto& bhe)
646 for (
auto const& in_out_component_id :
647 bhe.inflow_outflow_bc_component_ids)
649 if (bhe.use_python_bcs ||
650 this->_process_data._use_server_communication)
657 bcs.addBoundaryCondition(
659 get_global_bhe_bc_indices(
660 bhe.getBHEInflowDirichletBCNodesAndComponents(
661 bc_top_node_id, bc_bottom_node_id,
662 in_out_component_id.first)),
669 "The Python Boundary Condition was switched on, "
670 "but the data object does not exist! ");
682 get_global_bhe_bc_indices_with_bhe_idx(
684 {{{bc_top_node_id, in_out_component_id.first},
686 in_out_component_id.second}}}));
691 bcs.addBoundaryCondition(
693 get_global_bhe_bc_indices(
694 bhe.getBHEInflowDirichletBCNodesAndComponents(
695 bc_top_node_id, bc_bottom_node_id,
696 in_out_component_id.first)),
697 [&bhe](
double const T,
double const t) {
698 return bhe.updateFlowRateAndTemperature(T,
704 auto const bottom_nodes_and_components =
705 bhe.getBHEBottomDirichletBCNodesAndComponents(
707 in_out_component_id.first,
708 in_out_component_id.second);
710 if (bottom_nodes_and_components &&
716 bcs.addBoundaryCondition(
718 get_global_bhe_bc_indices(
719 {{{bc_bottom_node_id,
720 in_out_component_id.first},
722 in_out_component_id.second}}})));
724 else if (bottom_nodes_and_components &&
731 get_global_bhe_bc_indices_with_bhe_idx(
733 {{{bc_bottom_node_id, in_out_component_id.first},
735 in_out_component_id.second}}}));
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Global vector based on Eigen vector.
std::size_t searchByNodeIDs(const std::vector< std::size_t > &nodes)
Marks all elements connecting to any of the given nodes.
bool isAxiallySymmetric() const
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
std::size_t getID() const
Get id of the mesh.
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
std::size_t getNumberOfElements() const
Get the number of elements.
std::tuple< double, std::vector< double >, std::vector< double >, std::vector< int >, std::vector< double > > dataframe_network
bool isOverriddenTespy() const
virtual std::tuple< bool, bool, std::vector< double >, std::vector< double > > tespySolver(double, std::vector< double > const &, std::vector< double > const &) const
virtual void serverCommunicationPostTimestep(double, double, std::vector< double > const &, std::vector< double > const &, std::vector< double > const &) const
virtual std::tuple< std::vector< double >, std::vector< double > > serverCommunicationPreTimestep(double, double, std::vector< double > const &, std::vector< double > const &, std::vector< double > const &) const
bool isOverriddenServerCommunicationPostTimestep() const
bool isOverriddenServerCommunicationPreTimestep() const
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_soil_nodes
const BHEMeshData _bheMeshData
void assembleConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b) override
std::vector< std::unique_ptr< MeshLib::MeshSubset const > > _mesh_subset_BHE_nodes
std::vector< std::tuple< std::size_t, GlobalIndexType, GlobalIndexType > > _vec_bottom_BHE_node_indices
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id) override
std::vector< std::unique_ptr< HeatTransportBHELocalAssemblerInterface > > _local_assemblers
HeatTransportBHEProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const ¶meters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, HeatTransportBHEProcessData &&process_data, SecondaryVariableCollection &&secondary_variables, BHEMeshData &&bhe_mesh_data)
HeatTransportBHEProcessData _process_data
std::vector< std::size_t > _bhes_element_ids
NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &x) override
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
std::vector< std::tuple< std::size_t, GlobalIndexType, GlobalIndexType > > _vec_top_BHE_node_indices
AssembledMatrixCache _asm_mat_cache
void constructDofTable() override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void algebraicBcConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &xdot, int const process_id, GlobalMatrix &M, GlobalMatrix &K, GlobalVector &b)
void assembleWithJacobianConcreteProcess(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector &b, GlobalMatrix &Jac) override
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, const double t, const double dt, int const process_id) override
void createBHEBoundaryConditionTopBottom(std::vector< std::vector< MeshLib::Node * > > const &all_bhe_nodes)
std::vector< std::size_t > _soil_element_ids
virtual void computeSecondaryVariable(std::size_t const mesh_item_id, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_prev, int const process_id)
std::vector< BoundaryConditionCollection > _boundary_conditions
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
std::vector< std::size_t > const & getActiveElementIDs() const
VectorMatrixAssembler _global_assembler
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
MathLib::MatrixSpecifications getMatrixSpecifications(const int process_id) const override
Handles configuration of several secondary variables from the project file.
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, double const t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *M, GlobalMatrix *K, GlobalVector *b)
void assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalVector *b, GlobalMatrix *Jac)
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
void setMatrixSparsity(MATRIX &matrix, SPARSITY_PATTERN const &sparsity_pattern)
constexpr ranges::views::view_closure ids
For an element of a range view return its id.
PropertyVector< int > const * materialIDs(Mesh const &mesh)
@ BY_COMPONENT
Ordering data by component type.
std::unique_ptr< BHEInflowDirichletBoundaryCondition< BHEUpdateCallback > > createBHEInflowDirichletBoundaryCondition(std::pair< GlobalIndexType, GlobalIndexType > &&in_out_global_indices, BHEUpdateCallback bhe_update_callback)
std::unique_ptr< BHEBottomDirichletBoundaryCondition > createBHEBottomDirichletBoundaryCondition(std::pair< GlobalIndexType, GlobalIndexType > &&in_out_global_indices)
void createLocalAssemblers(std::vector< MeshLib::Element * > const &mesh_elements, NumLib::LocalToGlobalIndexMap const &dof_table, std::vector< std::unique_ptr< LocalAssemblerInterface > > &local_assemblers, NumLib::IntegrationOrder const integration_order, ExtraCtorArgs &&... extra_ctor_args)
std::unique_ptr< BHEInflowPythonBoundaryCondition< BHEType > > createBHEInflowPythonBoundaryCondition(std::pair< GlobalIndexType, GlobalIndexType > &&in_out_global_indices, BHEType &bhe, BHEInflowPythonBoundaryConditionPythonSideInterface &py_bc_object)
static void executeSelectedMemberOnDereferenced(Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
static void executeSelectedMemberDereferenced(Object &object, Method method, Container const &container, std::vector< std::size_t > const &active_container_ids, Args &&... args)
void assemble(const double t, double const dt, std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id, GlobalMatrix *const M, GlobalMatrix *const K, GlobalVector *const b, std::vector< NumLib::LocalToGlobalIndexMap const * > const &dof_tables, VectorMatrixAssembler &global_assembler, VectorOfLocalAssemblers const &local_assemblers, std::vector< std::size_t > const &active_element_ids)
const double _weighting_factor
const bool _use_algebraic_bc
std::vector< std::vector< MeshLib::Node * > > BHE_nodes
std::vector< std::vector< MeshLib::Element * > > BHE_elements
std::vector< int > BHE_mat_IDs
std::vector< bool > mass_lumping_soil_elements
MeshLib::PropertyVector< int > const * _mesh_prop_materialIDs
BHEInflowPythonBoundaryConditionPythonSideInterface * py_bc_object
Python object computing BC values.
std::unordered_map< int, int > _map_materialID_to_BHE_ID
std::vector< BHE::BHETypes > _vec_BHE_property
const bool _use_server_communication
AlgebraicBCSetting const _algebraic_BC_Setting