From 3099299d3fb62370f9a304aa6b9fb44d0d9e4e7b Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 13 Jul 2021 14:26:50 +0200
Subject: [PATCH] [python] Adapt FVProblem and add PorousMediumFlowProblem
---
dumux/python/CMakeLists.txt | 1 +
dumux/python/common/fvproblem.hh | 86 ++++++++-----
dumux/python/porousmediumflow/CMakeLists.txt | 3 +
dumux/python/porousmediumflow/problem.hh | 119 ++++++++++++++++++
python/dumux/CMakeLists.txt | 1 +
python/dumux/common/__init__.py | 16 +--
python/dumux/porousmediumflow/CMakeLists.txt | 3 +
python/dumux/porousmediumflow/__init__.py | 36 ++++++
test/python/test_explicit_transport_cctpfa.py | 5 +-
test/python/test_fvproblem.py | 4 +-
10 files changed, 235 insertions(+), 39 deletions(-)
create mode 100644 dumux/python/porousmediumflow/CMakeLists.txt
create mode 100644 dumux/python/porousmediumflow/problem.hh
create mode 100644 python/dumux/porousmediumflow/CMakeLists.txt
create mode 100644 python/dumux/porousmediumflow/__init__.py
diff --git a/dumux/python/CMakeLists.txt b/dumux/python/CMakeLists.txt
index ee2d01c436..4065a887e9 100644
--- a/dumux/python/CMakeLists.txt
+++ b/dumux/python/CMakeLists.txt
@@ -3,3 +3,4 @@ add_subdirectory(common)
add_subdirectory(discretization)
add_subdirectory(io)
add_subdirectory(material)
+add_subdirectory(porousmediumflow)
diff --git a/dumux/python/common/fvproblem.hh b/dumux/python/common/fvproblem.hh
index 5da3adc737..d4f9c8285c 100644
--- a/dumux/python/common/fvproblem.hh
+++ b/dumux/python/common/fvproblem.hh
@@ -42,12 +42,11 @@ namespace Dumux::Python {
* \ingroup Common
* \brief A C++ wrapper for a Python problem
*/
-template<class GridGeometry_, class PrimaryVariables, class SpatialParams_>
+template<class GridGeometry_, class PrimaryVariables, bool enableInternalDirichletConstraints_>
class FVProblem
{
public:
using GridGeometry = GridGeometry_;
- using SpatialParams = SpatialParams_;
using Scalar = typename PrimaryVariables::value_type;
using NumEqVector = Dune::FieldVector<Scalar, PrimaryVariables::dimension>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
@@ -61,25 +60,37 @@ public:
using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::dimension>;
FVProblem(std::shared_ptr<const GridGeometry> gridGeometry,
- std::shared_ptr<const SpatialParams> spatialParams,
pybind11::object pyProblem)
: gridGeometry_(gridGeometry)
- , spatialParams_(spatialParams)
, pyProblem_(pyProblem)
- {}
-
- std::string name() const
+ , name_("python_problem")
+ , paramGroup_("")
{
- return pyProblem_.attr("name").template cast<std::string>();
+ if (pybind11::hasattr(pyProblem_, "name"))
+ name_ = pyProblem.attr("name")().template cast<std::string>();
+
+ if (pybind11::hasattr(pyProblem_, "paramGroup"))
+ paramGroup_ = pyProblem.attr("paramGroup")().template cast<std::string>();
}
+ const std::string& name() const
+ { return name_; }
+
+ const std::string& paramGroup() const
+ { return paramGroup_; }
+
BoundaryTypes boundaryTypes(const Element &element,
const SubControlVolume &scv) const
{
if constexpr (!isBox)
DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scv) called for cell-centered method.");
else
- return pyProblem_.attr("boundaryTypes")(element, scv).template cast<BoundaryTypes>();
+ {
+ if (pybind11::hasattr(pyProblem_, "boundaryTypes"))
+ return pyProblem_.attr("boundaryTypes")(element, scv).template cast<BoundaryTypes>();
+ else
+ return pyProblem_.attr("boundaryTypesAtPos")(scv.dofPosition()).template cast<BoundaryTypes>();
+ }
}
BoundaryTypes boundaryTypes(const Element &element,
@@ -88,7 +99,12 @@ public:
if constexpr (isBox)
DUNE_THROW(Dune::InvalidStateException, "boundaryTypes(..., scvf) called for box method.");
else
- return pyProblem_.attr("boundaryTypes")(element, scvf).template cast<BoundaryTypes>();
+ {
+ if (pybind11::hasattr(pyProblem_, "boundaryTypes"))
+ return pyProblem_.attr("boundaryTypes")(element, scvf).template cast<BoundaryTypes>();
+ else
+ return pyProblem_.attr("boundaryTypesAtPos")(scvf.ipGlobal()).template cast<BoundaryTypes>();
+ }
}
PrimaryVariables dirichlet(const Element &element,
@@ -97,7 +113,12 @@ public:
if constexpr (!isBox)
DUNE_THROW(Dune::InvalidStateException, "dirichlet(scv) called for cell-centered method.");
else
- return pyProblem_.attr("dirichlet")(element, scv).template cast<PrimaryVariables>();
+ {
+ if (pybind11::hasattr(pyProblem_, "dirichlet"))
+ return pyProblem_.attr("dirichlet")(element, scv).template cast<PrimaryVariables>();
+ else
+ return pyProblem_.attr("dirichletAtPos")(scv.dofPosition()).template cast<PrimaryVariables>();
+ }
}
PrimaryVariables dirichlet(const Element &element,
@@ -106,7 +127,12 @@ public:
if constexpr (isBox)
DUNE_THROW(Dune::InvalidStateException, "dirichlet(scvf) called for box method.");
else
- return pyProblem_.attr("dirichlet")(element, scvf).template cast<PrimaryVariables>();
+ {
+ if (pybind11::hasattr(pyProblem_, "dirichlet"))
+ return pyProblem_.attr("dirichlet")(element, scvf).template cast<PrimaryVariables>();
+ else
+ return pyProblem_.attr("dirichletAtPos")(scvf.ipGlobal()).template cast<PrimaryVariables>();
+ }
}
template<class ElementVolumeVariables, class ElementFluxVariablesCache>
@@ -116,7 +142,10 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache,
const SubControlVolumeFace& scvf) const
{
- return pyProblem_.attr("neumann")(element, fvGeometry, scvf).template cast<NumEqVector>();
+ if (pybind11::hasattr(pyProblem_, "neumann"))
+ return pyProblem_.attr("neumann")(element, fvGeometry, scvf).template cast<NumEqVector>();
+ else
+ return pyProblem_.attr("neumannAtPos")(scvf.ipGlobal()).template cast<NumEqVector>();
}
template<class ElementVolumeVariables>
@@ -125,7 +154,10 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolume &scv) const
{
- return pyProblem_.attr("source")(element, fvGeometry, scv).template cast<NumEqVector>();
+ if (pybind11::hasattr(pyProblem_, "source"))
+ return pyProblem_.attr("source")(element, fvGeometry, scv).template cast<NumEqVector>();
+ else
+ return pyProblem_.attr("sourceAtPos")(scv.dofPosition()).template cast<NumEqVector>();
}
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
@@ -139,7 +171,10 @@ public:
const ElementVolumeVariables& elemVolVars,
const SubControlVolume& scv) const
{
- return pyProblem_.attr("scvPointSources")(element, fvGeometry, scv).template cast<NumEqVector>();
+ if (pybind11::hasattr(pyProblem_, "scvPointSources"))
+ return pyProblem_.attr("scvPointSources")(element, fvGeometry, scv).template cast<NumEqVector>();
+ else
+ return NumEqVector(0.0);
}
template<class Entity>
@@ -161,13 +196,8 @@ public:
return pyProblem_.attr("temperatureAtPos")(globalPos).template cast<Scalar>();
}
- const std::string paramGroup() const
- {
- return pyProblem_.attr("paramGroup")().template cast<std::string>();
- }
-
static constexpr bool enableInternalDirichletConstraints()
- { return false; } // TODO
+ { return enableInternalDirichletConstraints_; }
/*!
* \brief Add source term derivative to the Jacobian
@@ -180,20 +210,18 @@ public:
const VolumeVariables& volVars,
const SubControlVolume& scv) const
{
- pyProblem_.attr("addSourceDerivatives")(block, element, fvGeometry, scv).template cast<Scalar>();
+ if (pybind11::hasattr(pyProblem_, "addSourceDerivatives"))
+ pyProblem_.attr("addSourceDerivatives")(block, element, fvGeometry, scv);
}
-
const GridGeometry& gridGeometry() const
{ return *gridGeometry_; }
- const SpatialParams& spatialParams() const
- { return *spatialParams_; }
-
private:
std::shared_ptr<const GridGeometry> gridGeometry_;
- std::shared_ptr<const SpatialParams> spatialParams_;
pybind11::object pyProblem_;
+ std::string name_;
+ std::string paramGroup_;
};
// Python wrapper for the above FVProblem C++ class
@@ -204,11 +232,9 @@ void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options
using namespace Dune::Python;
using GridGeometry = typename Problem::GridGeometry;
- using SpatialParams = typename Problem::SpatialParams;
cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry,
- std::shared_ptr<const SpatialParams> spatialParams,
pybind11::object p){
- return std::make_shared<Problem>(gridGeometry, spatialParams, p);
+ return std::make_shared<Problem>(gridGeometry, p);
}));
cls.def_property_readonly("name", &Problem::name);
diff --git a/dumux/python/porousmediumflow/CMakeLists.txt b/dumux/python/porousmediumflow/CMakeLists.txt
new file mode 100644
index 0000000000..0b1c7f97be
--- /dev/null
+++ b/dumux/python/porousmediumflow/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_PYTHON_POROUSMEDIUMFLOW_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_PYTHON_POROUSMEDIUMFLOW_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/porousmediumflow)
diff --git a/dumux/python/porousmediumflow/problem.hh b/dumux/python/porousmediumflow/problem.hh
new file mode 100644
index 0000000000..30c72f3b40
--- /dev/null
+++ b/dumux/python/porousmediumflow/problem.hh
@@ -0,0 +1,119 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief TODO: docme!
+ */
+
+#ifndef DUMUX_PYTHON_POROUSMEDIUMFLOW_PROBLEM_HH
+#define DUMUX_PYTHON_POROUSMEDIUMFLOW_PROBLEM_HH
+
+
+#include <dune/python/pybind11/pybind11.h>
+
+#include <dumux/python/common/fvproblem.hh>
+
+namespace Dumux::Python {
+
+/*!
+ * \ingroup Common
+ * \brief A C++ wrapper for a Python PorousMediumFlow problem
+ */
+template<class GridGeometry_, class PrimaryVariables, class SpatialParams_,
+ bool enableInternalDirichletConstraints>
+class PorousMediumFlowProblem : public FVProblem<GridGeometry_, PrimaryVariables, enableInternalDirichletConstraints>
+{
+ using ParentType = FVProblem<GridGeometry_, PrimaryVariables, enableInternalDirichletConstraints>;
+public:
+ using GridGeometry = GridGeometry_;
+ using SpatialParams = SpatialParams_;
+ using Scalar = typename PrimaryVariables::value_type;
+ using NumEqVector = Dune::FieldVector<Scalar, PrimaryVariables::dimension>;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
+ static constexpr std::size_t numEq = static_cast<std::size_t>(PrimaryVariables::dimension);
+ using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::dimension>;
+
+ PorousMediumFlowProblem(std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<const SpatialParams> spatialParams,
+ pybind11::object pyProblem)
+ : ParentType(gridGeometry, pyProblem)
+ , spatialParams_(spatialParams)
+ {}
+
+ const SpatialParams& spatialParams() const
+ { return *spatialParams_; }
+
+private:
+ std::shared_ptr<const SpatialParams> spatialParams_;
+};
+
+// Python wrapper for the above FVProblem C++ class
+template<class Problem, class... options>
+void registerPorousMediumFlowProblem(pybind11::handle scope, pybind11::class_<Problem, options...> cls)
+{
+ using pybind11::operator""_a;
+ using namespace Dune::Python;
+
+ using GridGeometry = typename Problem::GridGeometry;
+ using SpatialParams = typename Problem::SpatialParams;
+ cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry,
+ std::shared_ptr<const SpatialParams> spatialParams,
+ pybind11::object p){
+ return std::make_shared<Problem>(gridGeometry, spatialParams, p);
+ }));
+
+ cls.def_property_readonly("name", &Problem::name);
+ cls.def_property_readonly("numEq", [](Problem&){ return Problem::numEq; });
+
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Vertex = typename GridView::template Codim<GridView::dimension>::Entity;
+
+ if constexpr (Problem::isBox)
+ {
+ using SCV = typename Problem::SubControlVolume;
+ cls.def("boundaryTypes", pybind11::overload_cast<const Element&, const SCV&>(&Problem::boundaryTypes, pybind11::const_), "element"_a, "scv"_a);
+ cls.def("dirichlet", pybind11::overload_cast<const Element&, const SCV&>(&Problem::dirichlet, pybind11::const_), "element"_a, "scv"_a);
+ }
+ else
+ {
+ using SCVF = typename Problem::SubControlVolumeFace;
+ cls.def("boundaryTypes", pybind11::overload_cast<const Element&, const SCVF&>(&Problem::boundaryTypes, pybind11::const_), "element"_a, "scvf"_a);
+ cls.def("dirichlet", pybind11::overload_cast<const Element&, const SCVF&>(&Problem::dirichlet, pybind11::const_), "element"_a, "scvf"_a);
+ }
+
+ cls.def("neumann", &Problem::template neumann<decltype(std::ignore), decltype(std::ignore)>);
+ cls.def("source", &Problem::template source<decltype(std::ignore)>);
+ cls.def("sourceAtPos", &Problem::sourceAtPos);
+ cls.def("initial", &Problem::template initial<Element>);
+ cls.def("initial", &Problem::template initial<Vertex>);
+ cls.def("extrusionFactor", &Problem::template extrusionFactor<decltype(std::ignore)>);
+ cls.def("gridGeometry", &Problem::gridGeometry);
+ cls.def("spatialParams", &Problem::spatialParams);
+}
+
+} // end namespace Dumux::Python
+
+#endif
diff --git a/python/dumux/CMakeLists.txt b/python/dumux/CMakeLists.txt
index 88b4fe39c8..db25915644 100644
--- a/python/dumux/CMakeLists.txt
+++ b/python/dumux/CMakeLists.txt
@@ -3,6 +3,7 @@ add_subdirectory(common)
add_subdirectory(discretization)
add_subdirectory(material)
add_subdirectory(io)
+add_subdirectory(porousmediumflow)
add_python_targets(dumux
__init__
diff --git a/python/dumux/common/__init__.py b/python/dumux/common/__init__.py
index 516a08206c..67a55de3a4 100644
--- a/python/dumux/common/__init__.py
+++ b/python/dumux/common/__init__.py
@@ -13,14 +13,15 @@ from dune.common.hashit import hashIt
# class MyProblem:
# ...
#
-def FVProblem(gridGeometry, spatialParams):
-
+def FVProblem(gridGeometry, enableInternalDirichletConstraints=False):
def createModule(numEq):
priVarType = "Dune::FieldVector<double, {}>".format(numEq)
ggType = gridGeometry._typeName
- spatialParamsType = spatialParams._typeName
- problemType = "Dumux::Python::FVProblem<{}, {}, {}>".format(ggType, priVarType, spatialParamsType)
- includes = gridGeometry._includes + spatialParams._includes + ["dumux/python/common/fvproblem.hh"]
+ enableIntDirConstraint = "true" if enableInternalDirichletConstraints else "false"
+ problemType = "Dumux::Python::FVProblem<{}, {}, {}>".format(
+ ggType, priVarType, enableIntDirConstraint
+ )
+ includes = gridGeometry._includes + ["dumux/python/common/fvproblem.hh"]
moduleName = "fvproblem_" + hashIt(problemType)
holderType = "std::shared_ptr<{}>".format(problemType)
generator = SimpleGenerator("FVProblem", "Dumux::Python")
@@ -31,7 +32,8 @@ def FVProblem(gridGeometry, spatialParams):
module = createModule(Cls.numEq)
def createFVProblem():
- return module.FVProblem(gridGeometry, spatialParams, Cls())
+ return module.FVProblem(gridGeometry, Cls())
+
return createFVProblem
return FVProblemDecorator
@@ -49,7 +51,7 @@ def BoundaryTypes(numEq=1):
moduleName = "boundarytypes_" + hashIt(typeName)
generator = SimpleGenerator("BoundaryTypes", "Dumux::Python")
module = generator.load(includes, typeName, moduleName)
- globals().update({cacheKey : module.BoundaryTypes})
+ globals().update({cacheKey: module.BoundaryTypes})
return globals()[cacheKey]()
diff --git a/python/dumux/porousmediumflow/CMakeLists.txt b/python/dumux/porousmediumflow/CMakeLists.txt
new file mode 100644
index 0000000000..9b6c72a388
--- /dev/null
+++ b/python/dumux/porousmediumflow/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_python_targets(porousmediumflow
+ __init__
+)
diff --git a/python/dumux/porousmediumflow/__init__.py b/python/dumux/porousmediumflow/__init__.py
new file mode 100644
index 0000000000..9db28887b5
--- /dev/null
+++ b/python/dumux/porousmediumflow/__init__.py
@@ -0,0 +1,36 @@
+from dune.generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+
+# A problem decorator generator for Python problems
+#
+# from dumux.common import PorousMediumFlowProblem
+# @PorousMediumFlowProblem(gridGeometry)
+# class MyProblem:
+# ...
+#
+def PorousMediumFlowProblem(gridGeometry, spatialParams, enableInternalDirichletConstraints=False):
+ def createModule(numEq):
+ priVarType = f"Dune::FieldVector<double, {numEq}>"
+ ggType = gridGeometry._typeName
+ spatialParamsType = spatialParams._typeName
+ enableIntDirConstraint = "true" if enableInternalDirichletConstraints else "false"
+ problemType = f"Dumux::Python::PorousMediumFlowProblem<{ggType}, {priVarType}, {spatialParamsType}, {enableIntDirConstraint}>"
+ includes = (
+ gridGeometry._includes
+ + spatialParams._includes
+ + ["dumux/python/porousmediumflow/problem.hh"]
+ )
+ moduleName = "fvproblem_" + hashIt(problemType)
+ generator = SimpleGenerator("PorousMediumFlowProblem", "Dumux::Python")
+ module = generator.load(includes, problemType, moduleName, holder="std::shared_ptr")
+ return module
+
+ def PorousMediumFlowProblemDecorator(Cls):
+ module = createModule(Cls.numEq)
+
+ def createPorousMediumFlowProblem():
+ return module.PorousMediumFlowProblem(gridGeometry, spatialParams, Cls())
+
+ return createPorousMediumFlowProblem
+
+ return PorousMediumFlowProblemDecorator
diff --git a/test/python/test_explicit_transport_cctpfa.py b/test/python/test_explicit_transport_cctpfa.py
index 06b18e5605..efa9282389 100755
--- a/test/python/test_explicit_transport_cctpfa.py
+++ b/test/python/test_explicit_transport_cctpfa.py
@@ -33,10 +33,13 @@ elementMapper = gridView.indexSet
# Define problem (inital/boundary condtions) #
##############################################
+
@FVProblem(gridGeometry)
class Problem:
numEq = 1
- name = "finitevolume"
+
+ def name(self):
+ return "finitevolume"
def boundaryTypes(self, element, scv):
bTypes = BoundaryTypes(self.numEq)
diff --git a/test/python/test_fvproblem.py b/test/python/test_fvproblem.py
index 74ecbf1057..cb3a64bb0a 100644
--- a/test/python/test_fvproblem.py
+++ b/test/python/test_fvproblem.py
@@ -29,7 +29,9 @@ gridGeometry.update()
@FVProblem(gridGeometry)
class Problem:
numEq = 2
- name = "python_problem"
+
+ def name(self):
+ return "python_problem"
def boundaryTypes(self, element, scv):
bTypes = BoundaryTypes(self.numEq)
--
GitLab