Merge pull request #2234 from TomFischer/PythonSourceTerms

Python source terms

Applications/CLI/CMakeLists.txt

@@ -39,7 +39,7 @@ if(OGS_USE_PYTHON)
     target_link_libraries(ogs_embedded_python PUBLIC pybind11::embed)
     target_compile_definitions(ogs_embedded_python PUBLIC OGS_USE_PYTHON)
     target_link_libraries(ogs_embedded_python PRIVATE
-        ProcessLibBoundaryConditionPythonModule)
+        ProcessLibBoundaryConditionPythonModule ProcessLibSourceTermPythonModule)
 
     target_link_libraries(ogs PRIVATE ogs_embedded_python)
 

Applications/CLI/ogs_embedded_python.cpp

@@ -13,12 +13,14 @@
 #include "ogs_embedded_python.h"
 
 #include "ProcessLib/BoundaryCondition/Python/PythonBoundaryConditionModule.h"
+#include "ProcessLib/SourceTerms/Python/PythonSourceTermModule.h"
 
 PYBIND11_EMBEDDED_MODULE(OpenGeoSys, m)
 {
     DBUG("Binding Python module OpenGeoSys.");
 
     ProcessLib::pythonBindBoundaryCondition(m);
+    ProcessLib::SourceTerms::Python::pythonBindSourceTerm(m);
 }
 
 #ifndef OGS_BUILD_SHARED_LIBS

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/c_Python.md

+Source term using a Python script to compute its values.
+
+See ProcessLib::SourceTerms::Python::PythonSourceTermPythonSideInterface for
+the Python-side source term interface.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_flush_stdout.md

+Determines if stdout should be flushed each time before and after passing
+control to the Python script. This might slow down operation, but ensures the
+right order of messages printed to stdout from Python and C++.

Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_source_term_object.md

+The object (variable name) in the provided Python script that will be used for
+this source term.

ProcessLib/CMakeLists.txt

@@ -32,6 +32,10 @@ if(OGS_USE_PYTHON)
     add_subdirectory(BoundaryCondition/Python)
     target_link_libraries(ProcessLib
         PUBLIC ProcessLibBoundaryConditionPython)
+
+    add_subdirectory(SourceTerms/Python)
+    target_link_libraries(ProcessLib
+        PUBLIC ProcessLibSourceTermPython)
 endif()
 
 if(OGS_INSITU)

ProcessLib/GroundwaterFlow/Tests.cmake

@@ -732,3 +732,26 @@ AddTest(
     python_laplace_eq_ref.vtu square_1e3_neumann_pcs_0_ts_1_t_1.000000.vtu pressure_expected pressure 4e-4 1e-16
 )
 
+AddTest(
+    NAME PythonSourceTermPoissonSinAXSinBYDirichlet_square_1e3
+    PATH Elliptic/square_1x1_GroundWaterFlow_Python
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e3_poisson_sin_x_sin_y.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS OGS_USE_PYTHON AND NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    DIFF_DATA
+    square_1x1_quad_1e3.vtu square_1e3_volumetricsourceterm_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure 0.7e-2 1e-16
+)
+
+AddTest(
+    NAME PythonSourceTermPoissonSinAXSinBYDirichlet_square_1e5
+    PATH Elliptic/square_1x1_GroundWaterFlow_Python
+    EXECUTABLE ogs
+    EXECUTABLE_ARGS square_1e5_poisson_sin_x_sin_y.prj
+    WRAPPER time
+    TESTER vtkdiff
+    REQUIREMENTS OGS_USE_PYTHON AND NOT (OGS_USE_LIS OR OGS_USE_MPI)
+    DIFF_DATA
+    square_1x1_quad_1e5.vtu square_1e5_volumetricsourceterm_pcs_0_ts_1_t_1.000000.vtu analytical_solution pressure 0.75e-4 1e-16
+)

ProcessLib/Process.cpp

@@ -187,9 +187,11 @@ void Process::assemble(const double t, GlobalVector const& x, GlobalMatrix& M,
 
     const auto pcs_id =
         (_coupled_solutions) != nullptr ? _coupled_solutions->process_id : 0;
+    // the last argument is for the jacobian, nullptr is for a unused jacobian
     _boundary_conditions[pcs_id].applyNaturalBC(t, x, K, b, nullptr);
 
-    _source_term_collections[pcs_id].integrate(t, b);
+    // the last argument is for the jacobian, nullptr is for a unused jacobian
+    _source_term_collections[pcs_id].integrate(t, x, b, nullptr);
 }
 
 void Process::assembleWithJacobian(const double t, GlobalVector const& x,

ProcessLib/SourceTerms/CreateSourceTerm.cpp

@@ -11,6 +11,9 @@
 
 #include "CreateNodalSourceTerm.h"
 #include "CreateVolumetricSourceTerm.h"
+#ifdef OGS_USE_PYTHON
+#include "Python/CreatePythonSourceTerm.h"
+#endif
 #include "SourceTerm.h"
 #include "SourceTermConfig.h"
 
@@ -41,6 +44,19 @@ std::unique_ptr<SourceTerm> createSourceTerm(
             *config.component_id);
     }
 
+    if (type == "Python")
+    {
+#ifdef OGS_USE_PYTHON
+        return ProcessLib::createPythonSourceTerm(
+            config.config, config.mesh, dof_table, mesh.getID(), variable_id,
+            *config.component_id, integration_order, shapefunction_order,
+            mesh.getDimension());
+#else
+        OGS_FATAL("OpenGeoSys has not been built with Python support.");
+#endif
+    }
+
+
     OGS_FATAL("Unknown source term type: `%s'.", type.c_str());
 }
 }  // namespace ProcessLib

ProcessLib/SourceTerms/NodalSourceTerm.cpp

@@ -37,7 +37,8 @@ NodalSourceTerm::NodalSourceTerm(
     }
 }
 
-void NodalSourceTerm::integrate(const double t, GlobalVector& b) const
+void NodalSourceTerm::integrate(const double t, GlobalVector const& /*x*/,
+                                GlobalVector& b, GlobalMatrix* /*jac*/) const
 {
     DBUG("Assemble NodalSourceTerm.");
 

ProcessLib/SourceTerms/NodalSourceTerm.h

@@ -22,7 +22,8 @@ public:
         const int variable_id, const int component_id,
         Parameter<double> const& parameter);
 
-    void integrate(const double t, GlobalVector& b) const override;
+    void integrate(const double t, GlobalVector const& x, GlobalVector& b,
+                   GlobalMatrix* jac) const override;
 
 private:
     std::size_t const _bulk_mesh_id;

ProcessLib/SourceTerms/Python/CMakeLists.txt

+add_library(ProcessLibSourceTermPython
+    CreatePythonSourceTerm.cpp
+    CreatePythonSourceTerm.h
+    PythonSourceTerm.cpp
+    PythonSourceTerm.h
+    PythonSourceTermLocalAssembler.h
+    PythonSourceTermPythonSideInterface.h)
+if(BUILD_SHARED_LIBS)
+    install(TARGETS ProcessLibSourceTermPython
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()
+
+target_compile_definitions(ProcessLibSourceTermPython
+    PUBLIC OGS_USE_PYTHON)
+
+target_link_libraries(ProcessLibSourceTermPython
+    PUBLIC BaseLib MathLib MeshLib NumLib logog
+    PRIVATE pybind11::pybind11)
+
+# For the embedded Python module
+add_library(ProcessLibSourceTermPythonModule
+    PythonSourceTermModule.cpp
+    PythonSourceTermModule.h)
+if(BUILD_SHARED_LIBS)
+    install(TARGETS ProcessLibSourceTermPythonModule
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR})
+endif()
+
+target_link_libraries(ProcessLibSourceTermPythonModule
+    PUBLIC
+    ProcessLibSourceTermPython
+    pybind11::pybind11)

ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "CreatePythonSourceTerm.h"
+
+#include <pybind11/pybind11.h>
+
+#include "BaseLib/ConfigTree.h"
+#include "MeshLib/Mesh.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/SourceTerms/SourceTerm.h"
+#include "PythonSourceTerm.h"
+
+namespace ProcessLib
+{
+std::unique_ptr<SourceTerm> createPythonSourceTerm(
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& source_term_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::size_t const bulk_mesh_id, int const variable_id,
+    int const component_id, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim)
+{
+    DBUG("Constructing PythonSourceTerm from config.");
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__type}
+    config.checkConfigParameter("type", "Python");
+
+    auto const source_term_object =
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__Python__source_term_object}
+        config.getConfigParameter<std::string>("source_term_object");
+    //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__Python__flush_stdout}
+    auto const flush_stdout = config.getConfigParameter("flush_stdout", false);
+
+    // Evaluate Python code in scope of main module
+    pybind11::object scope =
+        pybind11::module::import("__main__").attr("__dict__");
+
+    if (!scope.contains(source_term_object))
+        OGS_FATAL(
+            "Function `%s' is not defined in the python script file, or there "
+            "was no python script file specified.",
+            source_term_object.c_str());
+
+    auto* source_term = scope[source_term_object.c_str()]
+                            .cast<ProcessLib::SourceTerms::Python::
+                                      PythonSourceTermPythonSideInterface*>();
+
+    if (variable_id >= static_cast<int>(dof_table.getNumberOfVariables()) ||
+        component_id >= dof_table.getNumberOfVariableComponents(variable_id))
+    {
+        OGS_FATAL(
+            "Variable id or component id too high. Actual values: (%d, %d), "
+            "maximum values: (%d, %d).",
+            variable_id, component_id, dof_table.getNumberOfVariables(),
+            dof_table.getNumberOfVariableComponents(variable_id));
+    }
+
+    // In case of partitioned mesh the source_term could be empty, i.e. there is
+    // no source_term condition.
+#ifdef USE_PETSC
+    // This can be extracted to createSourceTerm() but then the config
+    // parameters are not read and will cause an error.
+    // TODO (naumov): Add a function to ConfigTree for skipping the tags of the
+    // subtree and move the code up in createSourceTerm().
+    if (source_term_mesh.getDimension() == 0 &&
+        source_term_mesh.getNumberOfNodes() == 0 &&
+        source_term_mesh.getNumberOfElements() == 0)
+    {
+        return nullptr;
+    }
+#endif  // USE_PETSC
+
+    return std::make_unique<ProcessLib::SourceTerms::Python::PythonSourceTerm>(
+        dof_table,
+        ProcessLib::SourceTerms::Python::PythonSourceTermData{
+            source_term, dof_table, bulk_mesh_id,
+            dof_table.getGlobalComponent(variable_id, component_id),
+            source_term_mesh},
+        integration_order, shapefunction_order, global_dim, flush_stdout);
+}
+
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MeshLib
+{
+class Mesh;
+}
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+}
+
+namespace ProcessLib
+{
+class SourceTerm;
+
+std::unique_ptr<SourceTerm> createPythonSourceTerm(
+    BaseLib::ConfigTree const& config, MeshLib::Mesh const& source_term_mesh,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
+    std::size_t const bulk_mesh_id, int const variable_id,
+    int const component_id, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim);
+
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "PythonSourceTerm.h"
+
+#include <pybind11/pybind11.h>
+#include <iostream>
+
+#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "PythonSourceTermLocalAssembler.h"
+
+namespace
+{
+//! Optionally flushes the standard output upon creation and destruction.
+//! Can be used to improve the debug output readability when printing debug
+//! messages both from OGS and from Python.
+class FlushStdoutGuard final
+{
+public:
+    //! Optionally flushes C++ stdout before running Python code.
+    explicit FlushStdoutGuard(bool const flush) : _flush(flush)
+    {
+        if (!flush)
+            return;
+
+        LOGOG_COUT << std::flush;
+    }
+
+    //! Optionally flushes Python's stdout after running Python code.
+    ~FlushStdoutGuard()
+    {
+        if (!_flush)
+            return;
+
+        using namespace pybind11::literals;
+        pybind11::print("end"_a = "", "flush"_a = true);
+    }
+
+private:
+    //! To flush or not to flush.
+    const bool _flush;
+};
+}  // anonymous namespace
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+PythonSourceTerm::PythonSourceTerm(
+    NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+    PythonSourceTermData&& source_term_data, unsigned const integration_order,
+    unsigned const shapefunction_order, unsigned const global_dim,
+    bool const flush_stdout)
+    : SourceTerm(source_term_dof_table),
+      _source_term_data(std::move(source_term_data)),
+      _flush_stdout(flush_stdout)
+{
+    std::vector<MeshLib::Node*> const& source_term_nodes =
+        _source_term_data.source_term_mesh.getNodes();
+    MeshLib::MeshSubset source_term_mesh_subset(
+        _source_term_data.source_term_mesh, source_term_nodes);
+
+    // Create local DOF table from the source term mesh subset for the given
+    // variable and component id.
+    _dof_table_source_term =
+        _source_term_data.dof_table_bulk.deriveBoundaryConstrainedMap(
+            std::move(source_term_mesh_subset));
+
+    createLocalAssemblers<PythonSourceTermLocalAssembler>(
+        global_dim, _source_term_data.source_term_mesh.getElements(),
+        *_dof_table_source_term, shapefunction_order, _local_assemblers,
+        _source_term_data.source_term_mesh.isAxiallySymmetric(),
+        integration_order, _source_term_data);
+}
+
+void PythonSourceTerm::integrate(const double t, const GlobalVector& x,
+                                 GlobalVector& b, GlobalMatrix* Jac) const
+{
+    FlushStdoutGuard guard(_flush_stdout);
+
+    GlobalExecutor::executeMemberOnDereferenced(
+        &PythonSourceTermLocalAssemblerInterface::assemble, _local_assemblers,
+        *_dof_table_source_term, t, x, b, Jac);
+}
+
+}  // namespace Python
+}  // namespace SourceTerms
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/PythonSourceTerm.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "ProcessLib/SourceTerms/SourceTerm.h"
+#include "PythonSourceTermLocalAssemblerInterface.h"
+#include "PythonSourceTermPythonSideInterface.h"
+
+namespace ProcessLib
+{
+class LocalToGlobalIndexMap;
+}
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Groups data used by source terms, in particular by the local assemblers.
+struct PythonSourceTermData final
+{
+    //! Python object computing source term values.
+    PythonSourceTermPythonSideInterface* source_term_object;
+
+    //! DOF table of the entire domain.
+    NumLib::LocalToGlobalIndexMap const& dof_table_bulk;
+
+    //! Mesh ID of the entire domain.
+    std::size_t const bulk_mesh_id;
+
+    //! Global component ID of the (variable, component) to which this source
+    //! term is applied.
+    int const global_component_id;
+
+    //! The source term mesh, i.e., the (sub-) domain of this source term.
+    const MeshLib::Mesh& source_term_mesh;
+};
+
+//! A source term whose values are computed by a Python script.
+class PythonSourceTerm final : public ProcessLib::SourceTerm
+{
+public:
+    explicit PythonSourceTerm(
+        NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+        PythonSourceTermData&& source_term_data,
+        unsigned const integration_order, unsigned const shapefunction_order,
+        unsigned const global_dim, bool const flush_stdout);
+
+    void integrate(const double t, GlobalVector const& x, GlobalVector& b,
+                   GlobalMatrix* jac) const override;
+
+private:
+    //! Auxiliary data.
+    PythonSourceTermData _source_term_data;
+
+    //! Local dof table for the boundary mesh.
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table_source_term;
+
+    //! Local assemblers for all elements of the source term mesh.
+    std::vector<std::unique_ptr<PythonSourceTermLocalAssemblerInterface>>
+        _local_assemblers;
+
+    //! Whether or not to flush standard output before and after each call to
+    //! Python code. Ensures right order of output messages and therefore
+    //! simplifies debugging.
+    bool const _flush_stdout;
+};
+
+}  // namespace Python
+}  // namespace SourceTerms
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "PythonSourceTerm.h"
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "PythonSourceTermLocalAssemblerInterface.h"
+#include "PythonSourceTermPythonSideInterface.h"
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Can be thrown to indicate that a member function is not overridden in a
+//! derived class (in particular, if a Python class inherits from a C++ class).
+struct MethodNotOverriddenInDerivedClassException
+{
+};
+
+template <typename NodalRowVectorType>
+struct IntegrationPointData final
+{
+    IntegrationPointData(NodalRowVectorType const& N_,
+                         double const& integration_weight_)
+        : N(N_), integration_weight(integration_weight_)
+    {
+    }
+
+    NodalRowVectorType const N;
+    double const integration_weight;
+
+    EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+          unsigned GlobalDim>
+class PythonSourceTermLocalAssembler final
+    : public PythonSourceTermLocalAssemblerInterface
+{
+public:
+    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+    using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+    using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+
+    explicit PythonSourceTermLocalAssembler(
+        MeshLib::Element const& e,
+        std::size_t const /*local_matrix_size*/,
+        bool is_axially_symmetric,
+        unsigned const integration_order,
+        PythonSourceTermData const& data)
+        : _data(data),
+          _element(e),
+          _is_axially_symmetric(is_axially_symmetric),
+          _integration_method(integration_order)
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        auto const shape_matrices =
+            initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                              IntegrationMethod, GlobalDim>(
+                _element, _is_axially_symmetric, _integration_method);
+
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
+        {
+            _ip_data.emplace_back(
+                shape_matrices[ip].N,
+                _integration_method.getWeightedPoint(ip).getWeight() *
+                    shape_matrices[ip].integralMeasure *
+                    shape_matrices[ip].detJ);
+        }
+    }
+
+    void assemble(std::size_t const source_term_element_id,
+                  NumLib::LocalToGlobalIndexMap const& dof_table_source_term,
+                  double const t, const GlobalVector& x, GlobalVector& b,
+                  GlobalMatrix* Jac) override
+    {
+        using ShapeMatricesType =
+            ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+        using FemType =
+            NumLib::TemplateIsoparametric<ShapeFunction, ShapeMatricesType>;
+        FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(
+            &_element));
+
+        unsigned const num_integration_points =
+            _integration_method.getNumberOfPoints();
+        auto const num_var = _data.dof_table_bulk.getNumberOfVariables();
+        auto const num_nodes = ShapeFunction::NPOINTS;
+        auto const num_comp_total =
+            _data.dof_table_bulk.getNumberOfComponents();
+
+        auto const& bulk_node_ids_map =
+            *_data.source_term_mesh.getProperties()
+                 .template getPropertyVector<std::size_t>("bulk_node_ids");
+
+        // gather primary variables
+        typename ShapeMatricesType::template MatrixType<ShapeFunction::NPOINTS,
+                                                       Eigen::Dynamic>
+            primary_variables_mat(num_nodes, num_comp_total);
+        for (int var = 0; var < num_var; ++var)
+        {
+            auto const num_comp =
+                _data.dof_table_bulk.getNumberOfVariableComponents(var);
+            for (int comp = 0; comp < num_comp; ++comp)
+            {
+                auto const global_component =
+                    dof_table_source_term.getGlobalComponent(var, comp);
+
+                for (unsigned element_node_id = 0; element_node_id < num_nodes;
+                     ++element_node_id)
+                {
+                    auto const* const node = _element.getNode(element_node_id);
+                    auto const boundary_node_id = node->getID();
+                    auto const bulk_node_id =
+                        bulk_node_ids_map[boundary_node_id];
+                    MeshLib::Location loc{_data.bulk_mesh_id,
+                                          MeshLib::MeshItemType::Node,
+                                          bulk_node_id};
+                    auto const dof_idx =
+                        _data.dof_table_bulk.getGlobalIndex(loc, var, comp);
+                    if (dof_idx == NumLib::MeshComponentMap::nop)
+                    {
+                        // TODO extend Python BC to mixed FEM ansatz functions
+                        OGS_FATAL(
+                            "No d.o.f. found for (node=%d, var=%d, comp=%d).  "
+                            "That might be due to the use of mixed FEM ansatz "
+                            "functions, which is currently not supported by "
+                            "the implementation of Python BCs. That excludes, "
+                            "e.g., the HM process.",
+                            bulk_node_id, var, comp);
+                    }
+                    primary_variables_mat(element_node_id, global_component) =
+                        x[dof_idx];
+                }
+            }
+        }
+
+        NodalRowVectorType local_rhs = Eigen::VectorXd::Zero(num_nodes);
+        NodalMatrixType local_Jac =
+            Eigen::MatrixXd::Zero(num_nodes, num_nodes * num_comp_total);
+        std::vector<double> prim_vars_data(num_comp_total);
+        auto prim_vars = MathLib::toVector(prim_vars_data);
+
+        for (unsigned ip = 0; ip < num_integration_points; ip++)
+        {
+            auto const& ip_data = _ip_data[ip];
+            auto const& N = ip_data.N;
+            auto const& coords = fe.interpolateCoordinates(N);
+            // Assumption: all primary variables have same shape functions.
+            prim_vars = N * primary_variables_mat;
+
+            auto const& w = ip_data.integration_weight;
+            auto const flux_dflux =
+                _data.source_term_object->getFlux(t, coords, prim_vars_data);
+            auto const& flux = flux_dflux.first;
+            auto const& dflux = flux_dflux.second;
+            local_rhs.noalias() += N * (flux * w);
+
+            if (static_cast<int>(dflux.size()) != num_comp_total)
+            {
+                // This strict check is technically mandatory only if a
+                // Jacobian is assembled. However, it is done as a
+                // consistency check also for cases without Jacobian
+                // assembly.
+                OGS_FATAL(
+                    "The Python source term must return the derivative of "
+                    "the flux w.r.t. each primary variable. %d components "
+                    "expected. %d components returned from Python.",
+                    num_comp_total, dflux.size());
+            }
+
+            if (Jac)
+            {
+                for (int comp = 0; comp < num_comp_total; ++comp)
+                {
+                    auto const top = 0;
+                    auto const left = comp * num_nodes;
+                    auto const width = num_nodes;
+                    auto const height = num_nodes;
+                    // The assignement -= takes into account the sign convention
+                    // of 1st-order in time ODE systems in OpenGeoSys.
+                    local_Jac.block(top, left, width, height).noalias() -=
+                        ip_data.N.transpose() * (dflux[comp] * w) * N;
+                }
+            }
+        }
+
+        auto const& indices_specific_component =
+            dof_table_source_term(source_term_element_id,
+                                  _data.global_component_id)
+                .rows;
+        b.add(indices_specific_component, local_rhs);
+
+        if (Jac)
+        {
+            // only assemble a block of the Jacobian, not the whole local matrix
+            auto const indices_all_components = NumLib::getIndices(
+                source_term_element_id, dof_table_source_term);
+            MathLib::RowColumnIndices<GlobalIndexType> rci{
+                indices_specific_component, indices_all_components};
+            Jac->add(rci, local_Jac);
+        }
+    }
+
+private:
+    PythonSourceTermData const& _data;
+    MeshLib::Element const& _element;
+    bool _is_axially_symmetric;
+
+    IntegrationMethod const _integration_method;
+    std::vector<
+        IntegrationPointData<NodalRowVectorType>,
+        Eigen::aligned_allocator<IntegrationPointData<NodalRowVectorType>>>
+        _ip_data;
+};
+
+}  // namespace Python
+}  // namespace SourceTerms
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h

+
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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
+{
+namespace SourceTerms
+{
+namespace Python
+{
+
+class PythonSourceTermLocalAssemblerInterface
+{
+public:
+    virtual void assemble(
+        std::size_t const source_term_element_id,
+        NumLib::LocalToGlobalIndexMap const& source_term_dof_table,
+        double const t, const GlobalVector& x, GlobalVector& b,
+        GlobalMatrix* Jac) = 0;
+
+    virtual ~PythonSourceTermLocalAssemblerInterface() = default;
+};
+
+}  // namespace Python
+}  // namespace SourceTerms
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "PythonSourceTermModule.h"
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include "PythonSourceTermPythonSideInterface.h"
+
+namespace ProcessLib
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Trampoline class allowing methods of class
+//! PythonSourceTermPythonSideInterface to be overridden on the Python
+//! side. Cf. https://pybind11.readthedocs.io/en/stable/advanced/classes.html
+class PythonSourceTermPythonSideInterfaceTrampoline
+    : public PythonSourceTermPythonSideInterface
+{
+public:
+    using PythonSourceTermPythonSideInterface::
+        PythonSourceTermPythonSideInterface;
+
+    std::pair<double, std::vector<double>> getFlux(
+        double t, std::array<double, 3> const& x,
+        std::vector<double> const& primary_variables) const override
+    {
+        using Ret = std::pair<double, std::vector<double>>;
+        PYBIND11_OVERLOAD_PURE(Ret, PythonSourceTermPythonSideInterface,
+                               getFlux, t, x, primary_variables);
+    }
+};
+
+void pythonBindSourceTerm(pybind11::module& m)
+{
+    namespace py = pybind11;
+
+    py::class_<PythonSourceTermPythonSideInterface,
+               PythonSourceTermPythonSideInterfaceTrampoline>
+        pybc(m, "SourceTerm");
+
+    pybc.def(py::init());
+
+    pybc.def("getFlux", &PythonSourceTermPythonSideInterface::getFlux);
+}
+
+}  // namespace Python
+}  // namespace SourceTerms
+}  // namespace ProcessLib

ProcessLib/SourceTerms/Python/PythonSourceTermModule.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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
+{
+namespace SourceTerms
+{
+namespace Python
+{
+//! Creates Python bindings for the Python source term class.
+OGS_EXPORT_SYMBOL void pythonBindSourceTerm(pybind11::module& m);
+}  // namespace Python
+}  // namespace SourceTerms
+}  // namespace ProcessLib