OGS
HeatTransportBHEProcess.cpp
Go to the documentation of this file.
1
12
13#include <cassert>
14
22
23namespace ProcessLib
24{
25namespace HeatTransportBHE
26{
28 std::string name,
29 MeshLib::Mesh& mesh,
30 std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
31 std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
32 unsigned const integration_order,
33 std::vector<std::vector<std::reference_wrapper<ProcessVariable>>>&&
34 process_variables,
35 HeatTransportBHEProcessData&& process_data,
36 SecondaryVariableCollection&& secondary_variables)
37 : Process(std::move(name), mesh, std::move(jacobian_assembler), parameters,
38 integration_order, std::move(process_variables),
39 std::move(secondary_variables)),
40 _process_data(std::move(process_data)),
41 _bheMeshData(getBHEDataInMesh(mesh))
42{
43 if (_bheMeshData.BHE_mat_IDs.size() !=
45 {
47 "The number of the given BHE properties ({:d}) are not consistent "
48 "with the number of BHE groups in the mesh ({:d}).",
51 }
52
53 auto material_ids = MeshLib::materialIDs(mesh);
54 if (material_ids == nullptr)
55 {
56 OGS_FATAL("Not able to get material IDs! ");
57 }
58
60
61 // create a map from a material ID to a BHE ID
62 for (int i = 0; i < static_cast<int>(_bheMeshData.BHE_mat_IDs.size()); i++)
63 {
64 // fill in the map structure
66 i;
67 }
68}
69
71{
72 // Create single component dof in every of the mesh's nodes.
74 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
75
76 //
77 // Soil temperature variable defined on the whole mesh.
78 //
80 std::make_unique<MeshLib::MeshSubset>(_mesh, _mesh.getNodes());
81 std::vector<MeshLib::MeshSubset> all_mesh_subsets{*_mesh_subset_soil_nodes};
82
83 std::vector<std::vector<MeshLib::Element*> const*> vec_var_elements;
84 vec_var_elements.push_back(&(_mesh.getElements()));
85
86 std::vector<int> vec_n_components{
87 1}; // one component for the soil temperature variable.
88
89 //
90 // BHE nodes with BHE type dependent number of variables.
91 //
92 int const n_BHEs = _process_data._vec_BHE_property.size();
93 assert(n_BHEs == static_cast<int>(_bheMeshData.BHE_mat_IDs.size()));
94 assert(n_BHEs == static_cast<int>(_bheMeshData.BHE_nodes.size()));
95 assert(n_BHEs == static_cast<int>(_bheMeshData.BHE_elements.size()));
96
97 // the BHE nodes need to be cherry-picked from the vector
98 for (int i = 0; i < n_BHEs; i++)
99 {
100 auto const number_of_unknowns =
101 visit([](auto const& bhe) { return bhe.number_of_unknowns; },
103 auto const& bhe_nodes = _bheMeshData.BHE_nodes[i];
104 auto const& bhe_elements = _bheMeshData.BHE_elements[i];
105
106 // All the BHE nodes have additional variables.
107 _mesh_subset_BHE_nodes.push_back(
108 std::make_unique<MeshLib::MeshSubset const>(_mesh, bhe_nodes));
109
110 std::generate_n(std::back_inserter(all_mesh_subsets),
111 // Here the number of components equals to the
112 // number of unknowns on the BHE
113 number_of_unknowns,
114 [&ms = _mesh_subset_BHE_nodes.back()]()
115 { return *ms; });
116
117 vec_n_components.push_back(number_of_unknowns);
118 vec_var_elements.push_back(&bhe_elements);
119 }
120
122 std::make_unique<NumLib::LocalToGlobalIndexMap>(
123 std::move(all_mesh_subsets),
124 vec_n_components,
125 vec_var_elements,
127
128 // in case of debugging the dof table, activate the following line
129 // std::cout << *_local_to_global_index_map << "\n";
130}
131
133 NumLib::LocalToGlobalIndexMap const& dof_table,
134 MeshLib::Mesh const& mesh,
135 unsigned const integration_order)
136{
137 // Quick access map to BHE's through element ids.
138 std::unordered_map<std::size_t, BHE::BHETypes*> element_to_bhe_map;
139 int const n_BHEs = _process_data._vec_BHE_property.size();
140 for (int i = 0; i < n_BHEs; i++)
141 {
142 auto const& bhe_elements = _bheMeshData.BHE_elements[i];
143 for (auto const& e : bhe_elements)
144 {
145 element_to_bhe_map[e->getID()] =
147 }
148 }
149
150 assert(mesh.getDimension() == 3);
153 mesh.getElements(), dof_table, _local_assemblers, element_to_bhe_map,
154 mesh.isAxiallySymmetric(), integration_order, _process_data);
155
156 // Create BHE boundary conditions for each of the BHEs
158}
159
161 const double t, double const dt, std::vector<GlobalVector*> const& x,
162 std::vector<GlobalVector*> const& xdot, int const process_id,
164{
165 DBUG("Assemble HeatTransportBHE process.");
166
167 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
168
169 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
170 dof_table = {std::ref(*_local_to_global_index_map)};
171 // Call global assembler for each local assembly item.
174 pv.getActiveElementIDs(), dof_table, t, dt, x, xdot, process_id, M, K,
175 b);
176}
177
179 const double t, double const dt, std::vector<GlobalVector*> const& x,
180 std::vector<GlobalVector*> const& xdot, int const process_id,
182{
183 DBUG("AssembleWithJacobian HeatTransportBHE process.");
184
185 std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
186 dof_table = {std::ref(*_local_to_global_index_map)};
187 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
188
189 // Call global assembler for each local assembly item.
192 _local_assemblers, pv.getActiveElementIDs(), dof_table, t, dt, x, xdot,
193 process_id, M, K, b, Jac);
194}
195
197 double const t, double const dt, std::vector<GlobalVector*> const& x,
198 GlobalVector const& x_dot, int const process_id)
199{
200 DBUG("Compute heat flux for HeatTransportBHE process.");
201
202 std::vector<NumLib::LocalToGlobalIndexMap const*> dof_tables;
203 dof_tables.reserve(x.size());
204 std::generate_n(std::back_inserter(dof_tables), x.size(),
205 [&]() { return _local_to_global_index_map.get(); });
206
207 ProcessLib::ProcessVariable const& pv = getProcessVariables(process_id)[0];
210 _local_assemblers, pv.getActiveElementIDs(), dof_tables, t, dt, x,
211 x_dot, process_id);
212}
213
214#ifdef OGS_USE_PYTHON
216 GlobalVector const& x)
217{
218 // if the process use python boundary condition
221
222 // Here the task is to get current time flowrate and flow temperature from
223 // TESPy and determine whether it converges.
224 auto const Tout_nodes_id =
226 const std::size_t n_bc_nodes = Tout_nodes_id.size();
227
228 for (std::size_t i = 0; i < n_bc_nodes; i++)
229 {
230 // read the T_out and store them in dataframe
232 x[Tout_nodes_id[i]];
233 }
234 // Transfer Tin and Tout to TESPy and return the results
235 auto const tespy_result = _process_data.py_bc_object->tespySolver(
237 std::get<1>(_process_data.py_bc_object->dataframe_network), // T_in
238 std::get<2>(_process_data.py_bc_object->dataframe_network)); // T_out
240 {
241 DBUG("Method `tespySolver' not overridden in Python script.");
242 }
243
244 // update the Tin and flow rate
245 for (std::size_t i = 0; i < n_bc_nodes; i++)
246 {
248 std::get<2>(tespy_result)[i];
250 std::get<3>(tespy_result)[i];
251 }
252 auto const tespy_has_converged = std::get<1>(tespy_result);
253 if (tespy_has_converged == true)
255
257}
258
260 std::vector<GlobalVector*> const& x, const double t, const double dt,
261 int const process_id)
262{
263 if (_process_data.py_bc_object == nullptr ||
265 {
266 return;
267 }
268
269 auto& [time, Tin_value, Tout_value, Tout_nodes_ids, flowrate] =
271
272 // We found the problem that time != t, but it always equals the last
273 // step. The following line is to correct this, although we do not use
274 // it for server communication.
275 time = t;
276
277 auto const& solution = *x[process_id];
278
279 // Iterate through each BHE
280 const std::size_t n_bc_nodes = Tout_nodes_ids.size();
281 for (std::size_t i = 0; i < n_bc_nodes; i++)
282 {
283 // read the T_out and store them in dataframe
284 Tout_value[i] = solution[Tout_nodes_ids[i]];
285 }
286
287 // Transfer T_out to server_Communication and get back T_in and flowrate
288 auto const server_communication_result =
290 Tout_value, flowrate);
293 {
294 DBUG("Method `serverCommunication' not overridden in Python script.");
295 }
296
297 auto const& [server_communication_Tin_value,
298 server_communication_flowrate] = server_communication_result;
299
300 std::copy(begin(server_communication_Tin_value),
301 end(server_communication_Tin_value),
302 begin(Tin_value));
303 std::copy(begin(server_communication_flowrate),
304 end(server_communication_flowrate),
305 begin(flowrate));
306}
307
309 std::vector<GlobalVector*> const& x, const double t, const double dt,
310 int const process_id)
311{
312 if (_process_data.py_bc_object == nullptr ||
314 {
315 return;
316 }
317
318 auto& [time, Tin_value, Tout_value, Tout_nodes_ids, flowrate] =
320
321 // We found the problem that time != t, but it always equals the last
322 // step. The following line is to correct this, although we do not use
323 // it for server communication.
324 time = t;
325
326 auto const& solution = *x[process_id];
327
328 // Iterate through each BHE
329 const std::size_t n_bc_nodes = Tout_nodes_ids.size();
330 for (std::size_t i = 0; i < n_bc_nodes; i++)
331 {
332 // read the T_out and store them in dataframe
333 Tout_value[i] = solution[Tout_nodes_ids[i]];
334 }
335
336 // Transfer T_out to server_Communication
338 Tout_value, flowrate);
341 {
342 DBUG("Method `serverCommunication' not overridden in Python script.");
343 }
344}
345#endif // OGS_USE_PYTHON
346
348 std::vector<std::vector<MeshLib::Node*>> const& all_bhe_nodes)
349{
350 const int process_id = 0;
351 auto& bcs = _boundary_conditions[process_id];
352
353 int const n_BHEs = static_cast<int>(_process_data._vec_BHE_property.size());
354
355 // for each BHE
356 for (int bhe_i = 0; bhe_i < n_BHEs; bhe_i++)
357 {
358 auto const& bhe_nodes = all_bhe_nodes[bhe_i];
359 // find the variable ID
360 // the soil temperature is 0-th variable
361 // the BHE temperature is therefore bhe_i + 1
362 const int variable_id = bhe_i + 1;
363
364 std::vector<MeshLib::Node*> bhe_boundary_nodes;
365
366 // cherry-pick the boundary nodes according to
367 // the number of connected line elements.
368 for (auto const& bhe_node : bhe_nodes)
369 {
370 // Count number of 1d elements connected with every BHE node.
371 auto const& connected_elements =
373 const std::size_t n_line_elements = std::count_if(
374 connected_elements.begin(), connected_elements.end(),
375 [](MeshLib::Element const* elem)
376 { return (elem->getDimension() == 1); });
377
378 if (n_line_elements == 1)
379 {
380 bhe_boundary_nodes.push_back(bhe_node);
381 }
382 }
383
384 if (bhe_boundary_nodes.size() != 2)
385 {
386 OGS_FATAL(
387 "Error!!! The BHE boundary nodes are not correctly found, "
388 "for every single BHE, there should be 2 boundary nodes.");
389 }
390
391 // For 1U, 2U, CXC, CXA type BHE, the node order in the boundary nodes
392 // vector should be rearranged according to its z coordinate in
393 // descending order. In these BHE types, the z coordinate on the top and
394 // bottom node is different. The BHE top node with a higher z coordinate
395 // should be placed at the first, while the BHE bottom node with a lower
396 // z coordinate should be placed at the second. For other horizontal BHE
397 // types e.g. 1P-type BHE, the z coordinate on the top and bottom node
398 // is identical. Thus the node order in the boundary nodes vector can
399 // not be rearranged according to its z coordinate. For these BHE types,
400 // the boundary node order is according to the default node id order in
401 // the model mesh.
402 // for 1P-type BHE
403 if ((*bhe_boundary_nodes[0])[2] == (*bhe_boundary_nodes[1])[2])
404 {
405 INFO(
406 "For 1P-type BHE, the BHE inflow and outflow "
407 "nodes are identified according to their mesh node id in "
408 "ascending order");
409 }
410 // for 1U, 2U, CXC, CXA type BHE
411 else
412 {
413 // swap the boundary nodes if the z coordinate of the
414 // first node is lower than it on the second node
415 if ((*bhe_boundary_nodes[0])[2] < (*bhe_boundary_nodes[1])[2])
416 {
417 std::swap(bhe_boundary_nodes[0], bhe_boundary_nodes[1]);
418 }
419 }
420
421 auto get_global_index =
422 [&](std::size_t const node_id, int const component)
423 {
424 return _local_to_global_index_map->getGlobalIndex(
426 variable_id, component);
427 };
428
429 auto get_global_bhe_bc_indices =
430 [&](std::array<
431 std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
432 nodes_and_components)
433 {
434 return std::make_pair(
435 get_global_index(nodes_and_components[0].first,
436 nodes_and_components[0].second),
437 get_global_index(nodes_and_components[1].first,
438 nodes_and_components[1].second));
439 };
440
441 auto createBCs =
442 [&, bc_top_node_id = bhe_boundary_nodes[0]->getID(),
443 bc_bottom_node_id = bhe_boundary_nodes[1]->getID()](auto& bhe)
444 {
445 for (auto const& in_out_component_id :
446 bhe.inflow_outflow_bc_component_ids)
447 {
448 if (bhe.use_python_bcs ||
450 // call BHEPythonBoundarycondition
451 {
452#ifdef OGS_USE_PYTHON
453 if (_process_data.py_bc_object) // the bc object exist
454 {
455 // apply the customized top, inflow BC.
456 bcs.addBoundaryCondition(
458 get_global_bhe_bc_indices(
459 bhe.getBHEInflowDirichletBCNodesAndComponents(
460 bc_top_node_id, bc_bottom_node_id,
461 in_out_component_id.first)),
462 bhe,
464 }
465 else
466 {
467 OGS_FATAL(
468 "The Python Boundary Condition was switched on, "
469 "but the data object does not exist! ");
470 }
471#else
472 OGS_FATAL(
473 "The Python Boundary Condition was switched off! "
474 "Not able to create Boundary Condition for BHE! ");
475#endif // OGS_USE_PYTHON
476 }
477 else
478 {
479 // Top, inflow, normal case
480 bcs.addBoundaryCondition(
482 get_global_bhe_bc_indices(
483 bhe.getBHEInflowDirichletBCNodesAndComponents(
484 bc_top_node_id, bc_bottom_node_id,
485 in_out_component_id.first)),
486 [&bhe](double const T, double const t) {
487 return bhe.updateFlowRateAndTemperature(T, t);
488 }));
489 }
490
491 auto const bottom_nodes_and_components =
492 bhe.getBHEBottomDirichletBCNodesAndComponents(
493 bc_bottom_node_id,
494 in_out_component_id.first,
495 in_out_component_id.second);
496
497 if (bottom_nodes_and_components)
498 {
499 // Bottom, outflow, all cases
500 bcs.addBoundaryCondition(
502 get_global_bhe_bc_indices(
503 {{{bc_bottom_node_id,
504 in_out_component_id.first},
505 {bc_bottom_node_id,
506 in_out_component_id.second}}})));
507 }
508 }
509 };
510 visit(createBCs, _process_data._vec_BHE_property[bhe_i]);
511 }
512}
513} // namespace HeatTransportBHE
514} // namespace ProcessLib
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:34
void DBUG(char const *fmt, Args const &... args)
Definition: Logging.h:29
Global vector based on Eigen vector.
Definition: EigenVector.h:28
bool isAxiallySymmetric() const
Definition: Mesh.h:131
std::vector< Node * > const & getNodes() const
Get the nodes-vector for the mesh.
Definition: Mesh.h:100
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:103
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:76
std::size_t getID() const
Get id of the mesh.
Definition: Mesh.h:115
std::vector< Element const * > const & getElementsConnectedToNode(std::size_t node_id) const
Definition: Mesh.cpp:238
std::tuple< double, std::vector< double >, std::vector< double >, std::vector< int >, std::vector< double > > dataframe_network
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
void computeSecondaryVariableConcrete(double const t, double const dt, std::vector< GlobalVector * > const &x, GlobalVector const &x_dot, int const process_id) override
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_soil_nodes
std::vector< std::unique_ptr< MeshLib::MeshSubset const > > _mesh_subset_BHE_nodes
void postTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
std::vector< std::unique_ptr< HeatTransportBHELocalAssemblerInterface > > _local_assemblers
NumLib::IterationResult postIterationConcreteProcess(GlobalVector const &x) override
HeatTransportBHEProcess(std::string name, MeshLib::Mesh &mesh, std::unique_ptr< ProcessLib::AbstractJacobianAssembler > &&jacobian_assembler, std::vector< std::unique_ptr< ParameterLib::ParameterBase > > const &parameters, unsigned const integration_order, std::vector< std::vector< std::reference_wrapper< ProcessVariable > > > &&process_variables, HeatTransportBHEProcessData &&process_data, SecondaryVariableCollection &&secondary_variables)
void preTimestepConcreteProcess(std::vector< GlobalVector * > const &x, const double t, const double dt, int const process_id) override
void assembleWithJacobianConcreteProcess(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, GlobalMatrix &Jac) override
void assembleConcreteProcess(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) override
void initializeConcreteProcess(NumLib::LocalToGlobalIndexMap const &dof_table, MeshLib::Mesh const &mesh, unsigned const integration_order) override
Process specific initialization called by initialize().
void createBHEBoundaryConditionTopBottom(std::vector< std::vector< MeshLib::Node * > > const &all_bhe_nodes)
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_dot, int const process_id)
std::vector< std::size_t > const & getActiveElementIDs() const
std::vector< BoundaryConditionCollection > _boundary_conditions
Definition: Process.h:367
std::unique_ptr< MeshLib::MeshSubset const > _mesh_subset_all_nodes
Definition: Process.h:329
MeshLib::Mesh & _mesh
Definition: Process.h:328
std::vector< std::reference_wrapper< ProcessVariable > > const & getProcessVariables(const int process_id) const
Definition: Process.h:146
VectorMatrixAssembler _global_assembler
Definition: Process.h:335
std::unique_ptr< NumLib::LocalToGlobalIndexMap > _local_to_global_index_map
Definition: Process.h:331
Handles configuration of several secondary variables from the project file.
void assemble(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, double const 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 assembleWithJacobian(std::size_t const mesh_item_id, LocalAssemblerInterface &local_assembler, std::vector< std::reference_wrapper< NumLib::LocalToGlobalIndexMap > > const &dof_tables, 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, GlobalMatrix &Jac)
IterationResult
Status flags telling the NonlinearSolver if an iteration succeeded.
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
static const double t
PropertyVector< int > const * materialIDs(Mesh const &mesh)
Definition: Mesh.cpp:264
@ 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, ExtraCtorArgs &&... extra_ctor_args)
BHEMeshData getBHEDataInMesh(MeshLib::Mesh const &mesh)
Definition: MeshUtils.cpp:51
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)
std::vector< std::vector< MeshLib::Node * > > BHE_nodes
Definition: MeshUtils.h:38
std::vector< std::vector< MeshLib::Element * > > BHE_elements
Definition: MeshUtils.h:37
BHEInflowPythonBoundaryConditionPythonSideInterface * py_bc_object
Python object computing BC values.