diff --git a/dumux/python/common/fvproblem.hh b/dumux/python/common/fvproblem.hh
index 2f9a1ec8bcdcc8b32292f9bdfe7f9c54672c26b6..5da3adc73700e7f94f4850b1bae381864163facc 100644
--- a/dumux/python/common/fvproblem.hh
+++ b/dumux/python/common/fvproblem.hh
@@ -42,11 +42,12 @@ namespace Dumux::Python {
* \ingroup Common
* \brief A C++ wrapper for a Python problem
*/
-template<class GridGeometry_, class PrimaryVariables>
+template<class GridGeometry_, class PrimaryVariables, class SpatialParams_>
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;
@@ -59,8 +60,12 @@ public:
static constexpr std::size_t numEq = static_cast<std::size_t>(PrimaryVariables::dimension);
using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::dimension>;
- FVProblem(std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object pyProblem)
- : gridGeometry_(gridGeometry), pyProblem_(pyProblem)
+ 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
@@ -128,6 +133,15 @@ public:
return pyProblem_.attr("sourceAtPos")(globalPos).template cast<NumEqVector>();
}
+ template<class ElementVolumeVariables>
+ NumEqVector scvPointSources(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume& scv) const
+ {
+ return pyProblem_.attr("scvPointSources")(element, fvGeometry, scv).template cast<NumEqVector>();
+ }
+
template<class Entity>
PrimaryVariables initial(const Entity& entity) const
{
@@ -142,11 +156,43 @@ public:
return pyProblem_.attr("extrusionFactor")(element, scv).template cast<Scalar>();
}
+ Scalar temperatureAtPos(const GlobalPosition& globalPos) const
+ {
+ 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
+
+ /*!
+ * \brief Add source term derivative to the Jacobian
+ * \note Only needed in case of analytic differentiation and solution dependent sources
+ */
+ template<class MatrixBlock, class VolumeVariables>
+ void addSourceDerivatives(MatrixBlock& block,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const VolumeVariables& volVars,
+ const SubControlVolume& scv) const
+ {
+ pyProblem_.attr("addSourceDerivatives")(block, element, fvGeometry, scv).template cast<Scalar>();
+ }
+
+
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_;
};
@@ -158,8 +204,11 @@ void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options
using namespace Dune::Python;
using GridGeometry = typename Problem::GridGeometry;
- cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object p){
- return std::make_shared<Problem>(gridGeometry, p);
+ 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);
diff --git a/dumux/python/material/spatialparams/spatialparams.hh b/dumux/python/material/spatialparams/spatialparams.hh
index e072d9bdd928b9c1e72136e832da4cae2a1d44d4..448a981491bcfca1747799614d79ee7be53d6b3b 100644
--- a/dumux/python/material/spatialparams/spatialparams.hh
+++ b/dumux/python/material/spatialparams/spatialparams.hh
@@ -24,22 +24,73 @@
#ifndef DUMUX_PYTHON_MATERIAL_SPATIAL_PARAMS_HH
#define DUMUX_PYTHON_MATERIAL_SPATIAL_PARAMS_HH
-
#include <dune/python/pybind11/pybind11.h>
#include <dune/python/pybind11/stl.h>
+#include <dumux/material/spatialparams/fv1p.hh>
+
namespace Dumux::Python {
+template<class GridGeometry, class Scalar, class PT>
+class FVSpatialParamsOneP
+: public Dumux::FVSpatialParamsOneP<GridGeometry, Scalar, FVSpatialParamsOneP<GridGeometry, Scalar, PT>>
+{
+ using ThisType = FVSpatialParamsOneP<GridGeometry, Scalar, PT>;
+ using ParentType = Dumux::FVSpatialParamsOneP<GridGeometry, Scalar, ThisType>;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using GlobalPosition = typename GridGeometry::GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+public:
+
+ using PermeabilityType = PT;
+
+ FVSpatialParamsOneP(std::shared_ptr<const GridGeometry> gridGeometry,
+ pybind11::object pySpatialParams)
+ : ParentType(gridGeometry)
+ , pySpatialParams_(pySpatialParams)
+ {}
+
+ template<class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ if (pybind11::hasattr(pySpatialParams_, "permeability"))
+ return pySpatialParams_.attr("permeability")(element, scv, elemSol).template cast<PermeabilityType>();
+ else
+ return pySpatialParams_.attr("permeabilityAtPos")(scv.center()).template cast<PermeabilityType>();
+ }
+
+ template<class ElementSolution>
+ Scalar porosity(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ if (pybind11::hasattr(pySpatialParams_, "porosity"))
+ return pySpatialParams_.attr("porosity")(element, scv, elemSol).template cast<Scalar>();
+ else
+ return pySpatialParams_.attr("porosityAtPos")(scv.center()).template cast<Scalar>();
+ }
+
+private:
+ pybind11::object pySpatialParams_;
+};
+
+
+
template <class SP, class... options>
-void registerSpatialParams(pybind11::handle scope,
- pybind11::class_<SP, options...> cls)
+void registerOnePSpatialParams(pybind11::handle scope,
+ pybind11::class_<SP, options...> cls)
{
using pybind11::operator""_a;
using GridGeometry = std::decay_t<decltype(std::declval<SP>().gridGeometry())>;
- cls.def(pybind11::init([](std::shared_ptr<GridGeometry> gridGeometry){
- return std::make_shared<SP>(gridGeometry);
- }), "gridGeometry"_a);
+ cls.def(pybind11::init([](std::shared_ptr<GridGeometry> gridGeometry, pybind11::object sp){
+ return std::make_shared<SP>(gridGeometry, sp);
+ }));
+
+ cls.def("gravity", &SP::gravity);
}
} // namespace Dumux::Python
diff --git a/python/dumux/common/__init__.py b/python/dumux/common/__init__.py
index 49c69da9f02af50c91ef77590885884130787fb9..516a08206cb06b45709008331541f1347f606d4c 100644
--- a/python/dumux/common/__init__.py
+++ b/python/dumux/common/__init__.py
@@ -13,13 +13,14 @@ from dune.common.hashit import hashIt
# class MyProblem:
# ...
#
-def FVProblem(gridGeometry):
+def FVProblem(gridGeometry, spatialParams):
def createModule(numEq):
priVarType = "Dune::FieldVector<double, {}>".format(numEq)
ggType = gridGeometry._typeName
- problemType = "Dumux::Python::FVProblem<{}, {}>".format(ggType, priVarType)
- includes = gridGeometry._includes + ["dumux/python/common/fvproblem.hh"]
+ spatialParamsType = spatialParams._typeName
+ problemType = "Dumux::Python::FVProblem<{}, {}, {}>".format(ggType, priVarType, spatialParamsType)
+ includes = gridGeometry._includes + spatialParams._includes + ["dumux/python/common/fvproblem.hh"]
moduleName = "fvproblem_" + hashIt(problemType)
holderType = "std::shared_ptr<{}>".format(problemType)
generator = SimpleGenerator("FVProblem", "Dumux::Python")
@@ -30,7 +31,7 @@ def FVProblem(gridGeometry):
module = createModule(Cls.numEq)
def createFVProblem():
- return module.FVProblem(gridGeometry, Cls())
+ return module.FVProblem(gridGeometry, spatialParams, Cls())
return createFVProblem
return FVProblemDecorator
diff --git a/python/dumux/material/__init__.py b/python/dumux/material/__init__.py
index dd20d3b3700effd3a40512aa63d34c7534b82248..48795f40e8c8ee5cc53e8e1f17aa44963148448e 100644
--- a/python/dumux/material/__init__.py
+++ b/python/dumux/material/__init__.py
@@ -1,3 +1,3 @@
from dumux.material.fluidsystems import FluidSystem
from dumux.material.components import Component
-from dumux.material.spatialparams import SpatialParams
+from dumux.material.spatialparams import OnePSpatialParams
diff --git a/python/dumux/material/spatialparams/__init__.py b/python/dumux/material/spatialparams/__init__.py
index 63472d5ca9f0cf4d5d2db113fb1cf3b4142eccf3..a2d3bbae36a982753cc70740fa371298cb635cef 100644
--- a/python/dumux/material/spatialparams/__init__.py
+++ b/python/dumux/material/spatialparams/__init__.py
@@ -1,14 +1,58 @@
from dune.generator.generator import SimpleGenerator
+from dumux.common import Property
+import numpy as np
-def SpatialParams(gridGeometry, scalar):
- moduleName = 'spatialparams'
- typeName = "Dumux::FVSpatialParamsOnePConstant<{}, {}>".format(gridGeometry._typeName, scalar._typeName)
- includes = ["dumux/material/spatialparams/fv1pconstant.hh"]
- includes += ["dumux/python/material/spatialparams/spatialparams.hh"]
- includes += gridGeometry._includes
- holderType = "std::shared_ptr<{}>".format(typeName)
- generator = SimpleGenerator("SpatialParams", "Dumux::Python")
- module = generator.load(includes, typeName, moduleName, options=[holderType])
+def OnePSpatialParams(*, gridGeometry, scalar=None):
+ moduleName = "spatialparams"
+ dim = gridGeometry.gridView.dimension
+ includes = gridGeometry._includes + [
+ "dumux/python/material/spatialparams/spatialparams.hh",
+ "dune/common/fmatrix.hh",
+ ]
- return module.SpatialParams(gridGeometry)
+ if scalar is None:
+ scalar = Property(type="double")
+
+ typeName = f"Dumux::Python::FVSpatialParamsOneP<{gridGeometry._typeName}, {scalar._typeName}, Dune::FieldMatrix<{scalar._typeName}, {dim}, {dim}>>"
+
+ def OnePSpatialParamsDecorator(Cls):
+ generator = SimpleGenerator("OnePSpatialParams", "Dumux::Python")
+ module = generator.load(includes, typeName, moduleName, holder="std::shared_ptr")
+
+ def maybeConvertScalarToMatrix(permeabilityValue):
+ if isinstance(permeabilityValue, float):
+ matrix = np.zeros(shape=(dim, dim))
+ np.fill_diagonal(matrix, permeabilityValue)
+ return matrix.tolist()
+ else:
+ return permeabilityValue
+
+ class Permeability:
+ def __init__(self, permeabilityFunction):
+ self.permeabilityFunction = permeabilityFunction
+
+ def __call__(self, element, scv, elemSol):
+ result = self.permeabilityFunction(element, scv, elemSol)
+ return maybeConvertScalarToMatrix(result)
+
+ class PermeabilityAtPos:
+ def __init__(self, permeabilityFunction):
+ self.permeabilityFunction = permeabilityFunction
+
+ def __call__(self, globalPos):
+ result = self.permeabilityFunction(globalPos)
+ return maybeConvertScalarToMatrix(result)
+
+ def createSpatialParams():
+ sp = Cls()
+ if hasattr(Cls, "permeability"):
+ Cls.permeability = Permeability(sp.permeability)
+ else:
+ Cls.permeabilityAtPos = PermeabilityAtPos(sp.permeabilityAtPos)
+ spatialParams = module.OnePSpatialParams(gridGeometry, sp)
+ return spatialParams
+
+ return createSpatialParams
+
+ return OnePSpatialParamsDecorator