diff --git a/Applications/CLI/CMakeLists.txt b/Applications/CLI/CMakeLists.txt
index f4e4a94213a15deeb51a0f63eb7f6b1a620a471d..74c05b7a45e9a407333d00bd035eefbc70145034 100644
--- a/Applications/CLI/CMakeLists.txt
+++ b/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)
diff --git a/Applications/CLI/ogs_embedded_python.cpp b/Applications/CLI/ogs_embedded_python.cpp
index 39b6e3e28182fe1fabe928e579a73dedd69f306b..b925f02e65e996c7e18a6c07e9d768473e9ae8b0 100644
--- a/Applications/CLI/ogs_embedded_python.cpp
+++ b/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
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/c_Python.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/c_Python.md
new file mode 100644
index 0000000000000000000000000000000000000000..dfd8a940617d636680e3ce482ecb203ca50d9c41
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/c_Python.md
@@ -0,0 +1,4 @@
+Source term using a Python script to compute its values.
+
+See ProcessLib::SourceTerms::Python::PythonSourceTermPythonSideInterface for
+the Python-side source term interface.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_flush_stdout.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_flush_stdout.md
new file mode 100644
index 0000000000000000000000000000000000000000..9ea47703dd005980f7a336584429e48635465a9d
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_flush_stdout.md
@@ -0,0 +1,3 @@
+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++.
diff --git a/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_source_term_object.md b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_source_term_object.md
new file mode 100644
index 0000000000000000000000000000000000000000..f36c8cd36ac7097d0fe9961ca75350142b0405dc
--- /dev/null
+++ b/Documentation/ProjectFile/prj/process_variables/process_variable/source_terms/source_term/Python/t_source_term_object.md
@@ -0,0 +1,2 @@
+The object (variable name) in the provided Python script that will be used for
+this source term.
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 61f2ee6746a5e216943913aa368515cdf45e6308..162ba4aacc5185ee79e8b4ad4f3c706984316932 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/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)
diff --git a/ProcessLib/GroundwaterFlow/Tests.cmake b/ProcessLib/GroundwaterFlow/Tests.cmake
index 13a983ac84126a5fa6a29df75ad98c5cb3189f15..da650d5a542909bec1c1dfe7ab5ab2f0ba8828e2 100644
--- a/ProcessLib/GroundwaterFlow/Tests.cmake
+++ b/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
+)
diff --git a/ProcessLib/Process.cpp b/ProcessLib/Process.cpp
index 152f2a3fba15e1f7014e318eca822955f75f29eb..a5e92bb5888ff4fad14b9f445e8ff73ab6b54fef 100644
--- a/ProcessLib/Process.cpp
+++ b/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,
diff --git a/ProcessLib/SourceTerms/CreateSourceTerm.cpp b/ProcessLib/SourceTerms/CreateSourceTerm.cpp
index 4e50e55f62343ea6eed0a69cb418af8ae9893279..b33cd640db6d86103e7bd2772411260ac2aca2d4 100644
--- a/ProcessLib/SourceTerms/CreateSourceTerm.cpp
+++ b/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
diff --git a/ProcessLib/SourceTerms/NodalSourceTerm.cpp b/ProcessLib/SourceTerms/NodalSourceTerm.cpp
index 04d8619f2ed413923c7fb21239bf2160de672cad..723b342672cfd1c04021b3591759a802c5d5608e 100644
--- a/ProcessLib/SourceTerms/NodalSourceTerm.cpp
+++ b/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.");
diff --git a/ProcessLib/SourceTerms/NodalSourceTerm.h b/ProcessLib/SourceTerms/NodalSourceTerm.h
index cc22ee10a8110f9f3db01bc2d60ddc538a965c90..4df7f58b07fdb9da6cdbb82bbcec93bdd973aa92 100644
--- a/ProcessLib/SourceTerms/NodalSourceTerm.h
+++ b/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;
diff --git a/ProcessLib/SourceTerms/Python/CMakeLists.txt b/ProcessLib/SourceTerms/Python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8ecddaea2e59a5715ef26188d8e27f4795415512
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/CMakeLists.txt
@@ -0,0 +1,32 @@
+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)
diff --git a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3707a3801e30d63d68c7193e29df0cbfbbd77cb4
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.cpp
@@ -0,0 +1,87 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..b7b2d6a4b7b328ba235a38ef7f6e12714b428f1a
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/CreatePythonSourceTerm.h
@@ -0,0 +1,38 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..733dd183140c9761507a12640b2220274651e9e0
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.cpp
@@ -0,0 +1,99 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTerm.h b/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ec77f3952f533ed64171d7a508f735d6256654c
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTerm.h
@@ -0,0 +1,79 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..9a8039ad6ce99ca501a78ba53069e4eb54b05da3
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
@@ -0,0 +1,234 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..f970bb39e4c799e06f3809705fa28fd7184d44ab
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssemblerInterface.h
@@ -0,0 +1,34 @@
+
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a556961c2d4a2c6882ae5a9d3a52d1e1f0c08093
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.cpp
@@ -0,0 +1,57 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermModule.h b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.h
new file mode 100644
index 0000000000000000000000000000000000000000..47b279fc0a34f28c3a400d44c488a6372d36b932
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermModule.h
@@ -0,0 +1,25 @@
+/**
+ * \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
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermPythonSideInterface.h b/ProcessLib/SourceTerms/Python/PythonSourceTermPythonSideInterface.h
new file mode 100644
index 0000000000000000000000000000000000000000..3314fb37d0a2aec4f5c9474205e691a6084af236
--- /dev/null
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermPythonSideInterface.h
@@ -0,0 +1,39 @@
+/**
+ * \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
+{
+//! Base class for source terms.
+//! This class will get Python bindings and is intended to be to be derived in
+//! Python.
+class PythonSourceTermPythonSideInterface
+{
+public:
+ /*!
+ * Computes the flux for the provided arguments (time, position of the node,
+ * primary variables at the node).
+ *
+ * \return flux Flux of the source term at that node and derivative of the
+ * flux w.r.t. all primary variables.
+ */
+ virtual std::pair<double, std::vector<double>> getFlux(
+ double /*t*/, std::array<double, 3> const& /*x*/,
+ std::vector<double> const& /*primary_variables*/) const = 0;
+
+ virtual ~PythonSourceTermPythonSideInterface() = default;
+};
+} // namespace Python
+} // namespace SourceTerms
+} // namespace ProcessLib
diff --git a/ProcessLib/SourceTerms/SourceTerm.h b/ProcessLib/SourceTerms/SourceTerm.h
index bbbd878a3e9a4e509a8954fc7909393293f20f30..de87c32c8431ef95ef70a44a6345fe5f86c1e9da 100644
--- a/ProcessLib/SourceTerms/SourceTerm.h
+++ b/ProcessLib/SourceTerms/SourceTerm.h
@@ -23,7 +23,8 @@ public:
{
}
- virtual void integrate(const double t, GlobalVector& b) const = 0;
+ virtual void integrate(const double t, GlobalVector const& x,
+ GlobalVector& b, GlobalMatrix* jac) const = 0;
virtual ~SourceTerm() = default;
diff --git a/ProcessLib/SourceTerms/SourceTermCollection.cpp b/ProcessLib/SourceTerms/SourceTermCollection.cpp
index 9636eae275a4460474d486082e03485bd099e18b..82bac27aeb295609eeee6ea49f52f8c365bec9eb 100644
--- a/ProcessLib/SourceTerms/SourceTermCollection.cpp
+++ b/ProcessLib/SourceTerms/SourceTermCollection.cpp
@@ -29,10 +29,11 @@ void SourceTermCollection::addSourceTermsForProcessVariables(
}
}
-void SourceTermCollection::integrate(const double t, GlobalVector& b) const
+void SourceTermCollection::integrate(const double t, GlobalVector const& x,
+ GlobalVector& b, GlobalMatrix* jac) const
{
for (auto const& st : _source_terms)
- st->integrate(t, b);
+ st->integrate(t, x, b, jac);
}
}
diff --git a/ProcessLib/SourceTerms/SourceTermCollection.h b/ProcessLib/SourceTerms/SourceTermCollection.h
index 77f8646c23abb1150f08956522373d991cd17e08..c1c26160aef2317699d751e40ee8abe778d628d4 100644
--- a/ProcessLib/SourceTerms/SourceTermCollection.h
+++ b/ProcessLib/SourceTerms/SourceTermCollection.h
@@ -23,7 +23,8 @@ public:
{
}
- void integrate(const double t, GlobalVector& b) const;
+ void integrate(const double t, GlobalVector const& x, GlobalVector& b,
+ GlobalMatrix* jac) const;
void addSourceTermsForProcessVariables(
std::vector<std::reference_wrapper<ProcessVariable>> const&
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
index 166405fa84ce45a1c0ab192e58406847675f7359..c05f41bdd3d342efbc027c66729ccd0d776ecda9 100644
--- a/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.cpp
@@ -47,8 +47,9 @@ VolumetricSourceTerm::VolumetricSourceTerm(
_volumetric_source_term);
}
-void VolumetricSourceTerm::integrate(const double t,
- GlobalVector& b) const
+void VolumetricSourceTerm::integrate(const double t, GlobalVector const& /*x*/,
+ GlobalVector& b,
+ GlobalMatrix* /*jac*/) const
{
DBUG("Assemble VolumetricSourceTerm.");
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTerm.h b/ProcessLib/SourceTerms/VolumetricSourceTerm.h
index 9f26f93d227613cccd89fa7f23632c8199050dcc..f3d785f310199d87d2cc086e72acaa07790089a7 100644
--- a/ProcessLib/SourceTerms/VolumetricSourceTerm.h
+++ b/ProcessLib/SourceTerms/VolumetricSourceTerm.h
@@ -27,7 +27,8 @@ public:
int const variable_id, int const component_id,
Parameter<double> const& volumetric_source_term);
- void integrate(const double t, GlobalVector& b) const;
+ void integrate(const double t, GlobalVector const& x, GlobalVector& b,
+ GlobalMatrix* jac) const;
private:
Parameter<double> const& _volumetric_source_term;
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/sin_x_sin_y_source_term.py b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/sin_x_sin_y_source_term.py
new file mode 100644
index 0000000000000000000000000000000000000000..60766042860fbb28f926e329e32f2743c299fb99
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/sin_x_sin_y_source_term.py
@@ -0,0 +1,24 @@
+
+import OpenGeoSys
+from math import pi, sin
+
+a = 2.0*pi
+b = 2.0*pi
+
+def solution(x, y):
+ return sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
+
+# - laplace(solution) = source term
+def laplace_solution(x, y):
+ return a*a * sin(a*x-pi/2.0) * sin(b*y-pi/2.0) + b*b * sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
+
+# source term for the benchmark
+class SinXSinYSourceTerm(OpenGeoSys.SourceTerm):
+ def getFlux(self, t, coords, primary_vars):
+ x, y, z = coords
+ value = laplace_solution(x,y)
+ Jac = [ 0.0 ]
+ return (value, Jac)
+
+# instantiate source term object referenced in OpenGeoSys' prj file
+sinx_siny_source_term = SinXSinYSourceTerm()
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_poisson_sin_x_sin_y.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_poisson_sin_x_sin_y.prj
new file mode 100644
index 0000000000000000000000000000000000000000..cf71bd3ad3a209e9949071b12fddd2dc42e4513e
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_poisson_sin_x_sin_y.prj
@@ -0,0 +1,132 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <meshes>
+ <mesh>square_1x1_quad_1e3.vtu</mesh>
+ <mesh>square_1x1_quad_1e3_geometry_ll.vtu</mesh>
+ <mesh>square_1x1_quad_1e3_geometry_lr.vtu</mesh>
+ <mesh>square_1x1_quad_1e3_geometry_ul.vtu</mesh>
+ <mesh>square_1x1_quad_1e3_geometry_ur.vtu</mesh>
+ <mesh>square_1x1_quad_1e3_entire_domain.vtu</mesh>
+ </meshes>
+ <python_script>sin_x_sin_y_source_term.py</python_script>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>square_1e3_volumetricsourceterm</prefix>
+ </output>
+ </time_loop>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>pressure_edge_points</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e3_geometry_ll</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e3_geometry_lr</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e3_geometry_ul</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e3_geometry_ur</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ <source_terms>
+ <source_term>
+ <mesh>square_1x1_quad_1e3_entire_domain</mesh>
+ <type>Python</type>
+ <source_term_object>sinx_siny_source_term</source_term_object>
+ </source_term>
+ </source_terms>
+ </process_variable>
+ </process_variables>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e5_poisson_sin_x_sin_y.prj b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e5_poisson_sin_x_sin_y.prj
new file mode 100644
index 0000000000000000000000000000000000000000..093c642323e109768c756eb2c65764150b268cb0
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e5_poisson_sin_x_sin_y.prj
@@ -0,0 +1,132 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <meshes>
+ <mesh>square_1x1_quad_1e5.vtu</mesh>
+ <mesh>square_1x1_quad_1e5_geometry_ll.vtu</mesh>
+ <mesh>square_1x1_quad_1e5_geometry_lr.vtu</mesh>
+ <mesh>square_1x1_quad_1e5_geometry_ul.vtu</mesh>
+ <mesh>square_1x1_quad_1e5_geometry_ur.vtu</mesh>
+ <mesh>square_1x1_quad_1e5_entire_domain.vtu</mesh>
+ </meshes>
+ <python_script>sin_x_sin_y_source_term.py</python_script>
+ <processes>
+ <process>
+ <name>GW23</name>
+ <type>GROUNDWATER_FLOW</type>
+ <integration_order>2</integration_order>
+ <hydraulic_conductivity>K</hydraulic_conductivity>
+ <process_variables>
+ <process_variable>pressure</process_variable>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable type="static" internal_name="darcy_velocity" output_name="v"/>
+ </secondary_variables>
+ </process>
+ </processes>
+ <time_loop>
+ <processes>
+ <process ref="GW23">
+ <nonlinear_solver>basic_picard</nonlinear_solver>
+ <convergence_criterion>
+ <type>DeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstol>1.e-6</abstol>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <output>
+ <variables>
+ <variable> pressure </variable>
+ <variable> v </variable>
+ </variables>
+ </output>
+ <time_stepping>
+ <type>SingleStep</type>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>square_1e5_volumetricsourceterm</prefix>
+ </output>
+ </time_loop>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>basic_picard</name>
+ <type>Picard</type>
+ <max_iter>10</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <lis>-i cg -p jacobi -tol 1e-16 -maxiter 10000</lis>
+ <eigen>
+ <solver_type>CG</solver_type>
+ <precon_type>DIAGONAL</precon_type>
+ <max_iteration_step>10000</max_iteration_step>
+ <error_tolerance>1e-16</error_tolerance>
+ </eigen>
+ <petsc>
+ <prefix>gw</prefix>
+ <parameters>-gw_ksp_type cg -gw_pc_type bjacobi -gw_ksp_rtol 1e-16 -gw_ksp_max_it 10000</parameters>
+ </petsc>
+ </linear_solver>
+ </linear_solvers>
+ <parameters>
+ <parameter>
+ <name>K</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0</value>
+ </parameter>
+ <parameter>
+ <name>pressure_edge_points</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ </parameters>
+ <process_variables>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e5_geometry_ll</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e5_geometry_lr</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e5_geometry_ul</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <mesh>square_1x1_quad_1e5_geometry_ur</mesh>
+ <type>Dirichlet</type>
+ <parameter>pressure_edge_points</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ <source_terms>
+ <source_term>
+ <mesh>square_1x1_quad_1e5_entire_domain</mesh>
+ <type>Python</type>
+ <source_term_object>sinx_siny_source_term</source_term_object>
+ </source_term>
+ </source_terms>
+ </process_variable>
+ </process_variables>
+</OpenGeoSysProject>
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
index 3e5bb2a80a8f6bbdc4b7776b696d7a6af0d43b51..c1cc186de8a9d70ba75dd3061a819ca6ca3dfe76 100644
--- a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1.gml
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:a3f52cd5a8058d967f8ead873e1afa323af39266e538bcd2c35af08683a81dad
-size 1083
+oid sha256:2f1dac68820a46835f6e89762fede4d9b6799bd125d18c624859639b95830808
+size 1123
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
index b56982101f6be9ed1cb17945723023e95d9b9f98..663c9c55ab81712cadb341c8c0ed292e06061245 100644
--- a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3.vtu
@@ -1,3 +1,3 @@
version https://git-lfs.github.com/spec/v1
-oid sha256:3653db864dd917e8d5be84d68ee2de13ba3485b1ef23de4e9894fb6749a1aa11
-size 25692
+oid sha256:d3ef47420d7071bfe1a8194048ab0de7b1ab53018814e08235a5c063c95d7609
+size 43824
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_entire_domain.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_entire_domain.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..982a4503f4a652f06d74371c1041b5deced18647
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_entire_domain.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5aa7d5132ac307d2bfb1bbe8ce737c6de309a3f494ad231f5eb512d7d4083d42
+size 155464
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ll.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ll.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2acfaf9fb2956770259901a765efd8d9d4938718
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ll.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:612e6cc7d6f49da22c69cf3b4987541f2c40c77e3abe14db8af4b20f8e966742
+size 1253
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_lr.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_lr.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..20373ccf7cbb79c1ceb4ff9ba5040dcdd0aea503
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_lr.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:83af233d08ebfbe5a8f987650f181161ae524fd7279f296d8bcd7aed131d9dd4
+size 1265
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ul.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ul.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..e45a14869950ad47325665d931c0f273fb083add
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ul.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9f09b957a59823c990360b50ed328f27390f8876a950987b75b31746ae748132
+size 1257
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ur.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ur.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..b79c7a7a1cffbf6cf565718ddcca94f52e11893c
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e3_geometry_ur.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:23183aa05576bac8c4e2a6255991c6ac69b4fc1fa32a7038ac60ff0dc7e2fb46
+size 1287
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a64c678a9b46fb204c2e9602e4df16d627dec94b
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ff1b3c1c89ed84a64f06112726f5653d0f5c74c2b646f58c8424bc5ba7d0ab99
+size 5678625
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_entire_domain.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_entire_domain.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..ca0d8d8fb8c30aef64adb7da27280589615c13d4
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_entire_domain.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:9c56618c5c8e4172801fead89bc46fdaa008be871992ac4427b213a24ef53fa3
+size 13489433
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ll.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ll.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..2acfaf9fb2956770259901a765efd8d9d4938718
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ll.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:612e6cc7d6f49da22c69cf3b4987541f2c40c77e3abe14db8af4b20f8e966742
+size 1253
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_lr.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_lr.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..3088e65f6a947ebc803304244610f234603f71f8
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_lr.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:a1646c78d38f7c6cbbed4b9a477baba2bcf254db55db52ae0d63188f77c5bd39
+size 1273
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ul.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ul.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..598e4dcf0c1a8698950d8d68e6bf1582405a6412
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ul.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:c1a75a0a8fe914c6baf5e2772825feee0924b414e319177c12b71c018c15abf2
+size 1261
diff --git a/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ur.vtu b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ur.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..df1cb1f6e84f916edc34872402e411c2373e8627
--- /dev/null
+++ b/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1x1_quad_1e5_geometry_ur.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:79d7684fd39a33fbe7a06f63060fb9dc67e6445938f3a6de5d6d5c70f685b307
+size 1295
diff --git a/web/content/docs/benchmarks/elliptic/poisson_equation.pandoc b/web/content/docs/benchmarks/elliptic/poisson_equation.pandoc
new file mode 100644
index 0000000000000000000000000000000000000000..0d920439850b933cc420168283f897931b071bbb
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/poisson_equation.pandoc
@@ -0,0 +1,194 @@
++++
+date = "2018-10-10T09:17:39+01:00"
+title = "Poisson equation using Python for source term specification"
+project = "Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_poisson_sin_x_sin_y.prj"
+author = "Tom Fischer"
+weight = 102
+
+aliases = [ "/docs/benchmarks/" ] # First benchmark page
+
+[menu]
+ [menu.benchmarks]
+ parent = "elliptic"
+
++++
+
+{{< data-link >}}
+
+## Poisson Equation
+
+The Poisson equation is:
+$$
+\begin{equation}
+- k\; \Delta p = Q \quad \text{in }\Omega
+\end{equation}$$
+w.r.t boundary conditions
+$$
+\eqalign{
+p(x) = g_D(x) &\quad \text{on }\Gamma_D,\cr
+k\;{\partial p(x) \over \partial n} = g_N(x) &\quad \text{on }\Gamma_N,
+}$$
+
+where $p$ could be the pressure, the subscripts $D$ and $N$ denote the
+Dirichlet- and Neumann-type boundary conditions, $n$ is the normal vector
+pointing outside of $\Omega$, and $\Gamma = \Gamma_D \cup \Gamma_N$ and
+$\Gamma_D \cap \Gamma_N = \emptyset$.
+
+## Problem specification and analytical solution
+
+We solve the Poisson equation on a square domain $\Omega = [0\times 1]^2$
+with $k = 1$ w.r.t. the specific boundary conditions:
+$$
+\eqalign{
+p(x,y) = 1 &\quad \text{on } (x=0,y=0) \subset \Gamma_D,\cr
+p(x,y) = 1 &\quad \text{on } (x=0,y=1) \subset \Gamma_D,\cr
+p(x,y) = 1 &\quad \text{on } (x=1,y=0) \subset \Gamma_D,\cr
+p(x,y) = 1 &\quad \text{on } (x=1,y=1) \subset \Gamma_D,\cr
+}$$
+and the source term is
+$$
+Q(x,y) = a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right)
++b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right)
+$$
+
+The analytical solution of (1) is
+$$
+p(x,y) = \sin\left(ax - \frac{\pi}{2}\right)
+ \sin\left(by - \frac{\pi}{2}\right).
+$$
+Since
+$$
+\frac{\partial^2 p}{\partial x} (x,y)
+ = - a^2 \sin\left(ax - \frac{\pi}{2}\right)
+ \sin\left(by - \frac{\pi}{2}\right)
+$$
+and
+$$
+\frac{\partial^2 p}{\partial y} (x,y)
+ = - b^2 \sin\left(ax - \frac{\pi}{2}\right)
+ \sin\left(by - \frac{\pi}{2}\right)
+$$
+it yields
+$$
+\Delta p(x,y)
+ = - a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right)
+ - b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right)
+$$
+and
+$$
+Q = - \Delta p(x,y)
+ = a^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right)
+ + b^2 \sin\left(ax - \frac{\pi}{2}\right) \sin\left(by - \frac{\pi}{2}\right).
+$$
+
+## Input files
+
+The project file is
+[`square_1e3_poisson_sin_x_sin_y.prj`](https://github.com/ufz/ogs/tree/master/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/square_1e3_poisson_sin_x_sin_y.prj). It describes the
+process to be solved and the related process variable together with their
+initial, boundary conditions and source terms. The definition of the source term
+$Q$ is in a [Python
+script](https://github.com/ufz/ogs/tree/master/Tests/Data/Elliptic/square_1x1_GroundWaterFlow_Python/sinx_siny_source_term.py).
+The script for setting the source terms is referenced in the project file as
+follows:
+```
+<python_script>sin_x_sin_y_source_term.py</python_script>
+```
+In the source term descripition
+```
+<source_term>
+ <mesh>square_1x1_quad_1e3_entire_domain</mesh>
+ <type>Python</type>
+ <source_term_object>sinx_siny_source_term</source_term_object>
+</source_term>
+```
+the domain is specified by the mesh-tag. The function $Q$ is defined by the
+Python object `sinx_sinx_source_term` that is created in the last line of the
+Python script:
+```python
+import OpenGeoSys
+from math import pi, sin
+
+a = 2.0*pi
+b = 2.0*pi
+
+def solution(x, y):
+ return - sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
+
+# - laplace(solution) = source term
+def laplace_solution(x, y):
+ return a*a * sin(a*x-pi/2.0) * sin(b*y-pi/2.0) + b*b * sin(a*x-pi/2.0) * sin(b*y-pi/2.0)
+
+# source term for the benchmark
+class SinXSinYSourceTerm(OpenGeoSys.SourceTerm):
+ def getFlux(self, t, coords, primary_vars):
+ x, y, z = coords
+ value = laplace_solution(x,y)
+ Jac = [ 0.0 ]
+ return (value, Jac)
+
+# instantiate source term object referenced in OpenGeoSys' prj file
+sinx_siny_source_term = SinXSinYSourceTerm()
+```
+
+## Running simulation
+
+To start the simulation (after successful compilation) run:
+```bash
+$ ogs square_1e3_poisson_sin_x_sin_y.prj
+```
+
+It will produce some output and write the computed result into the
+`square_1e3_volumetricsourceterm_pcs_0_ts_1_t_1.000000.vtu`, which can be
+directly visualized and analysed in paraview for example.
+
+The output on the console will be similar to the following on:
+```
+info: This is OpenGeoSys-6 version 6.1.0-1132-g00a6062a2.
+info: OGS started on 2018-10-10 09:22:17+0200.
+
+info: ConstantParameter: K
+info: ConstantParameter: p0
+info: ConstantParameter: pressure_edge_points
+info: Initialize processes.
+createPythonSourceTerm
+info: Solve processes.
+info: [time] Output of timestep 0 took 0.000695229 s.
+info: === Time stepping at step #1 and time 1 with step size 1
+info: [time] Assembly took 0.0100119 s.
+info: [time] Applying Dirichlet BCs took 0.000133991 s.
+info: ------------------------------------------------------------------
+info: *** Eigen solver computation
+info: -> solve with CG (precon DIAGONAL)
+info: iteration: 81/10000
+info: residual: 5.974447e-17
+
+info: ------------------------------------------------------------------
+info: [time] Linear solver took 0.00145817 s.
+info: [time] Iteration #1 took 0.0116439 s.
+info: [time] Solving process #0 took 0.011662 s in time step #1
+info: [time] Time step #1 took 0.0116858 s.
+info: [time] Output of timestep 1 took 0.000671864 s.
+info: The whole computation of the time stepping took 1 steps, in which
+ the accepted steps are 1, and the rejected steps are 0.
+
+info: [time] Execution took 0.0370049 s.
+info: OGS terminated on 2018-10-10 09:22:17+020
+```
+
+## Results and evaluation
+
+### Comparison of the numerical and analytical solutions
+
+{{< img src="../square_1e3_poisson_sin_x_sin_y_sourceterm_Pressure_PythonSourceTerm.png" >}}
+The above picture shows the numerical result. The solution conforms in the edges
+to the prescribed boundary conditions.
+{{< img src="../square_1e3_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png" >}}
+Since a coarse mesh ($32 \times 32$ elements) is used for the simulation the
+difference between the numerical and the analytical solution is relatively large.
+
+#### Comparison for higher resolution mesh ($316 \times 316$ elements)
+{{< img src="../square_1e5_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png" >}}
+The difference between the numerical and the analytical solution is much smaller
+than in the coarse mesh case.
+
diff --git a/web/content/docs/benchmarks/elliptic/square_1e3_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png b/web/content/docs/benchmarks/elliptic/square_1e3_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png
new file mode 100644
index 0000000000000000000000000000000000000000..9c4e18ab60faec0c22dea3b93e51412ec116858c
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/square_1e3_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:776901f0e2b1eac0114c682265bb95578a37fcc66616f040e3be97528c6682ac
+size 164143
diff --git a/web/content/docs/benchmarks/elliptic/square_1e3_poisson_sin_x_sin_y_sourceterm_Pressure_PythonSourceTerm.png b/web/content/docs/benchmarks/elliptic/square_1e3_poisson_sin_x_sin_y_sourceterm_Pressure_PythonSourceTerm.png
new file mode 100644
index 0000000000000000000000000000000000000000..fc2efbc0e2aed28239e5b16cab83ccfe7a1a02c3
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/square_1e3_poisson_sin_x_sin_y_sourceterm_Pressure_PythonSourceTerm.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:8f0535683ab8662903d1d8dcf38c3bf290d0925591a0ce5b2b99dfe3a0048567
+size 183103
diff --git a/web/content/docs/benchmarks/elliptic/square_1e5_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png b/web/content/docs/benchmarks/elliptic/square_1e5_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png
new file mode 100644
index 0000000000000000000000000000000000000000..dcfc2d318d1be110dc27a0a5780944685eb11c9f
--- /dev/null
+++ b/web/content/docs/benchmarks/elliptic/square_1e5_poisson_sin_x_sin_y_sourceterm_Diff_Pressure_AnalyticalSolution_PythonSourceTerm.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:5bede41f6571cb9d91e6fd619c2eb08757c06f772842ad1426e4c9cc6345de8c
+size 163087