Add BHEInflowPythonBoundaryCondition

Add BHEinflowPythonBoundaryCondition into following places
1) PythonBoundaryCondition
2) Pybind11_Embedded_Module
3) The CMakeList
4) Switch on the compiling flag OGS_USE_PYTHON when pipe network is present
Add OGS_USE_PYTHON option in HeatTransportBHE process cmakelist

Applications/CLI/ogs_embedded_python.cpp

@@ -7,11 +7,13 @@
  *
  */
 
+#include "ogs_embedded_python.h"
+
 #include <pybind11/embed.h>
-#include <logog/include/logog.hpp>
 
-#include "ogs_embedded_python.h"
+#include <logog/include/logog.hpp>
 
+#include "ProcessLib/BoundaryCondition/Python/BHEInflowPythonBoundaryConditionModule.h"
 #include "ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h"
 #include "ProcessLib/SourceTerms/Python/PythonSourceTermModule.h"
 
@@ -20,6 +22,7 @@ PYBIND11_EMBEDDED_MODULE(OpenGeoSys, m)
     DBUG("Binding Python module OpenGeoSys.");
 
     ProcessLib::pythonBindBoundaryCondition(m);
+    ProcessLib::bheInflowpythonBindBoundaryCondition(m);
     ProcessLib::SourceTerms::Python::pythonBindSourceTerm(m);
 }
 

ProcessLib/BoundaryCondition/Python/BHEInflowPythonBoundaryCondition.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+#ifdef OGS_USE_PYTHON
+#include <pybind11/pybind11.h>
+#endif  // OGS_USE_PYTHON
+#include <algorithm>
+#include <iostream>
+#include <logog/include/logog.hpp>
+#include <vector>
+
+#include "BHEInflowPythonBoundaryConditionPythonSideInterface.h"
+#include "BaseLib/Error.h"
+#include "NumLib/IndexValueVector.h"
+#include "ProcessLib/BoundaryCondition/BoundaryCondition.h"
+#include "ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h"
+#include "ProcessLib/HeatTransportBHE/BHE/BHETypes.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "PythonBoundaryConditionLocalAssembler.h"
+
+namespace ProcessLib
+{
+//! A boundary condition whose values are computed by a Python script.
+template <typename BHEType>
+class BHEInflowPythonBoundaryCondition final : public BoundaryCondition
+{
+public:
+    BHEInflowPythonBoundaryCondition(
+        std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
+        BHEType& bhe,
+        BHEInflowPythonBoundaryConditionPythonSideInterface& py_bc_object)
+        : _in_out_global_indices(std::move(in_out_global_indices)),
+          _bhe(bhe),
+          _py_bc_object(py_bc_object)
+    {
+        const auto g_idx_T_out = in_out_global_indices.second;
+
+        // store the bc node ids to BHE network dataframe
+        std::get<4>(_py_bc_object.dataframe_network).emplace_back(g_idx_T_out);
+    }
+
+    void getEssentialBCValues(
+        const double t, const GlobalVector& x,
+        NumLib::IndexValueVector<GlobalIndexType>& bc_values) const override
+    {
+        bc_values.ids.resize(1);
+        bc_values.values.resize(1);
+        auto const& data_exchange = _py_bc_object.dataframe_network;
+        // get the number of all boundary nodes
+        const std::size_t n_bc_nodes = std::get<4>(data_exchange).size();
+
+        // get T_in bc_id
+        bc_values.ids[0] = _in_out_global_indices.first;
+
+        // get T_out bc_id
+        auto const boundary_node_id = _in_out_global_indices.second;
+
+        // return T_in from currently BHE dataframe column 2,
+        // update flowrate and HeatTransferCoefficients for each BHE
+        for (std::size_t i = 0; i < n_bc_nodes; i++)
+        {
+            // auto pair_flag_value =
+            // _bc_data.bc_object->getDirichletBCValue(boundary_node_id);
+            auto const dataframe_node_id = std::get<4>(data_exchange);
+            auto const dataframe_Tin_val = std::get<2>(data_exchange);
+            auto const dataframe_BHE_flowrate = std::get<5>(data_exchange);
+            if (dataframe_node_id[i] == boundary_node_id)
+            {
+                bc_values.values[0] = dataframe_Tin_val[i];
+                _bhe.updateHeatTransferCoefficients(dataframe_BHE_flowrate[i]);
+                break;
+            }
+        }
+
+        // store the current time to network dataframe
+        std::get<1>(_py_bc_object.dataframe_network) = t;
+    }
+
+private:
+    std::pair<GlobalIndexType, GlobalIndexType> const _in_out_global_indices;
+    BHEType& _bhe;
+    BHEInflowPythonBoundaryConditionPythonSideInterface& _py_bc_object;
+};
+
+template <typename BHEType>
+std::unique_ptr<BHEInflowPythonBoundaryCondition<BHEType>>
+createBHEInflowPythonBoundaryCondition(
+    std::pair<GlobalIndexType, GlobalIndexType>&& in_out_global_indices,
+    BHEType& bhe,
+    BHEInflowPythonBoundaryConditionPythonSideInterface& py_bc_object)
+
+{
+    DBUG("Constructing BHEInflowPythonBoundaryCondition.");
+
+    // In case of partitioned mesh the boundary could be empty, i.e. there is no
+    // boundary condition.
+#ifdef USE_PETSC
+    // For this special boundary condition the boundary condition is not empty
+    // if the global indices are non-negative.
+    if (in_out_global_indices.first < 0 && in_out_global_indices.second < 0)
+    {
+        return nullptr;
+    }
+    // If only one of the global indices (in or out) is negative the
+    // implementation is not valid.
+    if (in_out_global_indices.first < 0 || in_out_global_indices.second < 0)
+    {
+        OGS_FATAL(
+            "The partition cuts the BHE into two independent parts. This "
+            "behaviour is not implemented.");
+    }
+#endif  // USE_PETSC
+    return std::make_unique<BHEInflowPythonBoundaryCondition<BHEType>>(
+        std::move(in_out_global_indices), bhe, py_bc_object);
+}
+}  // namespace ProcessLib

ProcessLib/BoundaryCondition/Python/BHEInflowPythonBoundaryConditionModule.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BHEInflowPythonBoundaryConditionModule.h"
+
+#include <pybind11/stl.h>
+
+#include "BHEInflowPythonBoundaryConditionPythonSideInterface.h"
+
+namespace ProcessLib
+{
+//! Trampoline class allowing methods of class
+//! PythonBoundaryConditionPythonSideInterface to be overridden on the Python
+//! side. Cf. https://pybind11.readthedocs.io/en/stable/advanced/classes.html
+class BHEInflowPythonBoundaryConditionPythonSideInterfaceTrampoline
+    : public BHEInflowPythonBoundaryConditionPythonSideInterface
+{
+public:
+    using BHEInflowPythonBoundaryConditionPythonSideInterface::
+        BHEInflowPythonBoundaryConditionPythonSideInterface;
+
+    std::tuple<bool, double, std::vector<double>, std::vector<double>,
+               std::vector<int>, std::vector<double>>
+    initializeDataContainer() const override
+    {
+        using Ret =
+            std::tuple<bool, double, std::vector<double>, std::vector<double>,
+                       std::vector<int>, std::vector<double>>;
+        PYBIND11_OVERLOAD(Ret,
+                          BHEInflowPythonBoundaryConditionPythonSideInterface,
+                          initializeDataContainer);
+    }
+
+    std::tuple<bool, bool, std::vector<double>> tespyThermalSolver(
+        double t,
+        std::vector<double> const& Tin_val,
+        std::vector<double> const& Tout_val) const override
+    {
+        using Ret = std::tuple<bool, bool, std::vector<double>>;
+        PYBIND11_OVERLOAD(Ret,
+                          BHEInflowPythonBoundaryConditionPythonSideInterface,
+                          tespyThermalSolver, t, Tin_val, Tout_val);
+    }
+
+    std::tuple<bool, std::vector<double>> tespyHydroSolver(
+        double t) const override
+    {
+        using Ret = std::tuple<bool, std::vector<double>>;
+        PYBIND11_OVERLOAD(Ret,
+                          BHEInflowPythonBoundaryConditionPythonSideInterface,
+                          tespyHydroSolver, t);
+    }
+};
+
+void bheInflowpythonBindBoundaryCondition(pybind11::module& m)
+{
+    namespace py = pybind11;
+
+    py::class_<BHEInflowPythonBoundaryConditionPythonSideInterface,
+               BHEInflowPythonBoundaryConditionPythonSideInterfaceTrampoline>
+        pybc(m, "BHENetwork");
+
+    pybc.def(py::init());
+
+    pybc.def("initializeDataContainer",
+             &BHEInflowPythonBoundaryConditionPythonSideInterface::
+                 initializeDataContainer);
+    pybc.def("tespyThermalSolver",
+             &BHEInflowPythonBoundaryConditionPythonSideInterface::
+                 tespyThermalSolver);
+    pybc.def(
+        "tespyHydroSolver",
+        &BHEInflowPythonBoundaryConditionPythonSideInterface::tespyHydroSolver);
+}
+
+}  // namespace ProcessLib

ProcessLib/BoundaryCondition/Python/BHEInflowPythonBoundaryConditionModule.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <pybind11/pybind11.h>
+
+#include "BaseLib/ExportSymbol.h"
+
+namespace ProcessLib
+{
+//! Creates BHE Inflow Python bindings for the Python BC class.
+OGS_EXPORT_SYMBOL void bheInflowpythonBindBoundaryCondition(
+    pybind11::module& m);
+}  // namespace ProcessLib

ProcessLib/BoundaryCondition/Python/BHEInflowPythonBoundaryConditionPythonSideInterface.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2019, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace ProcessLib
+{
+//! Base class for BHENetwork.
+//! This class will get Python bindings and is intended to be to be derived in
+//! Python.
+class BHEInflowPythonBoundaryConditionPythonSideInterface
+{
+public:
+    /*!
+     * Initialize network dataframe
+     * return a tuple (time, BHE inflow temperature, BHE outflow
+     * temperature, BHE outflow bc node id, BHE flowrate)
+     * set at that position and the parameters of the BHE network.
+     */
+    virtual std::tuple<bool,
+                       double /*time*/,
+                       std::vector<double> /*Tin_val*/,
+                       std::vector<double> /*Tout_val*/,
+                       std::vector<int> /*bc_out_ids*/,
+                       std::vector<double> /*BHE_flowrate*/>
+    initializeDataContainer() const
+    {
+        _overridden_essential = false;
+        return std::tuple<bool,
+                          double,
+                          std::vector<double>,
+                          std::vector<double>,
+                          std::vector<int>,
+                          std::vector<double>>{
+            false, std::numeric_limits<double>::quiet_NaN(), {}, {}, {}, {}};
+    }
+
+    /*!
+     * transfer BHE network dataframe to TESPy and get Tin from TESPy
+     *
+     * \return a tuple (time, BHE Tin and Tout value from TESPy)
+     * indicating if a natural BC shall be set at that position and the new
+     * inflow temperature of all BHEs
+     */
+    virtual std::tuple<bool, bool, std::vector<double>> tespyThermalSolver(
+        double /*t*/,
+        std::vector<double> const& /*Tin_val*/,
+        std::vector<double> const& /*Tout_val*/) const
+    {
+        _overridden_natural = false;
+        return std::tuple<bool, bool, std::vector<double>>{false, false, {}};
+    }
+
+    /*!
+     * call Tespy hydraulic solver to get flow velocity in each pipe
+     *
+     * \return a tuple (is_natural,  f_velocity ) indicating if a
+     * natural BC shall be set at that position and the flow velocity in each
+     * pipe of all BHEs
+     */
+    virtual std::tuple<bool, std::vector<double>> tespyHydroSolver(
+        double /*t*/) const
+    {
+        _overridden_natural = false;
+        return std::tuple<bool, std::vector<double>>{false, {}};
+    }
+
+    //! Tells if initializeDataContainer() has been overridden in the derived
+    //! class in Python.
+    //!
+    //! \pre initializeDataContainer() must already have been called
+    //! once.
+    bool isOverriddenEssential() const { return _overridden_essential; }
+
+    //! Tells if tespySolver() has been overridden in the derived class in
+    //! Python.
+    //!
+    //! \pre tespySolver() must already have been called once.
+    bool isOverriddenNatural() const { return _overridden_natural; }
+
+    // BHE network dataframe container
+    std::tuple<bool,
+               double,
+               std::vector<double>,
+               std::vector<double>,
+               std::vector<int>,
+               std::vector<double>>
+        dataframe_network;
+
+    virtual ~BHEInflowPythonBoundaryConditionPythonSideInterface() = default;
+
+private:
+    //! Tells if initializeDataContainer() has been overridden in the derived
+    //! class in Python.
+    mutable bool _overridden_essential = true;
+    //! Tells if tespySolver() has been overridden in the derived class in
+    //! Python.
+    mutable bool _overridden_natural = true;
+};
+}  // namespace ProcessLib

ProcessLib/BoundaryCondition/Python/CMakeLists.txt

 add_library(ProcessLibBoundaryConditionPython
             PythonBoundaryCondition.cpp PythonBoundaryCondition.h
             PythonBoundaryConditionLocalAssembler.h
-            PythonBoundaryConditionPythonSideInterface.h)
+            PythonBoundaryConditionPythonSideInterface.h
+            BHEInflowPythonBoundaryCondition.h
+            BHEInflowPythonBoundaryConditionPythonSideInterface.h)
 if(BUILD_SHARED_LIBS)
     install(TARGETS ProcessLibBoundaryConditionPython
             LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
@@ -21,7 +23,8 @@ target_link_libraries(ProcessLibBoundaryConditionPython
 
 # For the embedded Python module
 add_library(ProcessLibBoundaryConditionPythonModule
-            PythonBoundaryConditionModule.cpp PythonBoundaryConditionModule.h)
+            PythonBoundaryConditionModule.cpp PythonBoundaryConditionModule.h
+            BHEInflowPythonBoundaryConditionModule.cpp BHEInflowPythonBoundaryConditionModule.h)
 if(BUILD_SHARED_LIBS)
     install(TARGETS ProcessLibBoundaryConditionPythonModule
             LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})

ProcessLib/HeatTransportBHE/CMakeLists.txt

@@ -4,6 +4,13 @@ append_source_files(SOURCES BoundaryConditions)
 append_source_files(SOURCES LocalAssemblers)
 
 add_library(HeatTransportBHE ${SOURCES})
+
+if(OGS_USE_PYTHON)
+    target_link_libraries(HeatTransportBHE PUBLIC ProcessLib PRIVATE pybind11::pybind11)
+else()
+    target_link_libraries(HeatTransportBHE PUBLIC ProcessLib)
+endif()
+
 if(BUILD_SHARED_LIBS)
     install(TARGETS HeatTransportBHE
             LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})