Commit 790bdaee authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/spatialparams_python' into 'master'

Spatial parameters for Python bindings

Closes #1145

See merge request !3069
parents 6a70b718 97d568bb
Pipeline #15923 passed with stage
in 1 minute and 14 seconds
......@@ -35,6 +35,7 @@
#include <dumux/common/boundarytypes.hh>
#include <dumux/discretization/method.hh>
#include <dumux/python/common/boundarytypes.hh>
#include <dumux/python/common/fvspatialparams.hh>
namespace Dumux::Python {
......@@ -42,12 +43,13 @@ namespace Dumux::Python {
* \ingroup Common
* \brief A C++ wrapper for a Python problem
*/
template<class GridGeometry_, class PrimaryVariables, bool enableInternalDirichletConstraints_>
template<class GridGeometry_, class SpatialParams_, class PrimaryVariables, bool enableInternalDirichletConstraints_>
class FVProblem
{
public:
using GridGeometry = GridGeometry_;
using Scalar = typename PrimaryVariables::value_type;
using SpatialParams = SpatialParams_;
using Scalar = typename GridGeometry::GridView::ctype;
using NumEqVector = Dune::FieldVector<Scalar, PrimaryVariables::dimension>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GridGeometry::LocalView;
......@@ -60,11 +62,13 @@ 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)
, pyProblem_(pyProblem)
, name_("python_problem")
, paramGroup_("")
, spatialParams_(spatialParams)
{
if (pybind11::hasattr(pyProblem_, "name"))
name_ = pyProblem.attr("name")().template cast<std::string>();
......@@ -73,6 +77,11 @@ public:
paramGroup_ = pyProblem.attr("paramGroup")().template cast<std::string>();
}
FVProblem(std::shared_ptr<const GridGeometry> gridGeometry,
pybind11::object pyProblem)
: FVProblem(gridGeometry, std::make_shared<SpatialParams>(gridGeometry), pyProblem)
{}
const std::string& name() const
{ return name_; }
......@@ -183,19 +192,6 @@ public:
return pyProblem_.attr("initial")(entity).template cast<PrimaryVariables>();
}
template<class ElementSolution>
Scalar extrusionFactor(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
return pyProblem_.attr("extrusionFactor")(element, scv).template cast<Scalar>();
}
Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
return pyProblem_.attr("temperatureAtPos")(globalPos).template cast<Scalar>();
}
static constexpr bool enableInternalDirichletConstraints()
{ return enableInternalDirichletConstraints_; }
......@@ -214,14 +210,39 @@ public:
pyProblem_.attr("addSourceDerivatives")(block, element, fvGeometry, scv);
}
[[deprecated("extrusionFactorAtPos() should now be defined in the spatial params. This interface will be removed after 3.5.")]]
Scalar extrusionFactorAtPos(const GlobalPosition &globalPos, double defaultValue = 1.0) const
{ return 1.0; }
template<class ElementSolution>
[[deprecated("extrusionFactor() should now be defined in the spatial params. This interface will be removed after 3.5.")]]
Scalar extrusionFactor(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol,
double defaultValue = 1.0) const
{ return this->extrusionFactorAtPos(scv.center()); }
[[deprecated("temperatureAtPos() should now be defined in the spatial params. This interface will be removed after 3.5.")]]
Scalar temperatureAtPos(const GlobalPosition &globalPos, int defaultValue = 1) const
{ return 293.0; }
[[deprecated("temperature() should now be defined in the spatial params. This interface will be removed after 3.5.")]]
Scalar temperature(int defaultValue = 1) const
{ return this->temperature(GlobalPosition(0.0)); }
const GridGeometry& gridGeometry() const
{ return *gridGeometry_; }
//! Return a reference to the underlying spatial parameters
const SpatialParams& spatialParams() const
{ return *spatialParams_; }
private:
std::shared_ptr<const GridGeometry> gridGeometry_;
pybind11::object pyProblem_;
std::string name_;
std::string paramGroup_;
std::shared_ptr<const SpatialParams> spatialParams_;
};
// Python wrapper for the above FVProblem C++ class
......@@ -232,6 +253,12 @@ 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);
}));
cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry,
pybind11::object p){
return std::make_shared<Problem>(gridGeometry, p);
......@@ -263,7 +290,6 @@ void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options
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)>);
}
} // end namespace Dumux::Python
......
// -*- 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
* \ingroup Common
* \ingroup SpatialParameters
* \brief Basic spatial parameters to be used with finite-volume schemes.
*/
#ifndef DUMUX_PYTHON_COMMON_FVSPATIALPARAMS_HH
#define DUMUX_PYTHON_COMMON_FVSPATIALPARAMS_HH
#include <memory>
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dune/python/pybind11/pybind11.h>
namespace Dumux::Python {
/*!
* \ingroup Common
* \ingroup SpatialParameters
* \brief The base class for spatial parameters used with finite-volume schemes.
*/
template<class GridGeometry_>
class FVSpatialParams
{
public:
using GridGeometry = GridGeometry_;
using GridView = typename GridGeometry::GridView;
using Scalar = typename GridView::ctype;
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolume = typename GridGeometry::SubControlVolume;
static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethods::box;
static constexpr int dimWorld = GridView::dimensionworld;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using GravityVector = Dune::FieldVector<Scalar, dimWorld>;
FVSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry,
pybind11::object pySpatialParameters)
: gridGeometry_(gridGeometry)
, pySpatialParameters_(pySpatialParameters)
, gravity_(0.0)
{
if (getParam<bool>("Problem.EnableGravity"))
gravity_[dimWorld-1] = -9.81;
}
FVSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
: gridGeometry_(gridGeometry)
, pySpatialParameters_{}
, gravity_(0.0)
{
if (getParam<bool>("Problem.EnableGravity"))
gravity_[dimWorld-1] = -9.81;
}
/*!
* \brief Return how much the domain is extruded at a given sub-control volume.
*
* This means the factor by which a lower-dimensional (1D or 2D)
* entity needs to be expanded to get a full dimensional cell. The
* default is 1.0 which means that 1D problems are actually
* thought as pipes with a cross section of 1 m^2 and 2D problems
* are assumed to extend 1 m to the back.
*/
template<class ElementSolution>
Scalar extrusionFactor(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
if (pySpatialParameters_)
{
if (pybind11::hasattr(pySpatialParameters_, "extrusionFactor"))
return pySpatialParameters_.attr("extrusionFactor")(element, scv, elemSol).template cast<Scalar>();
else if (pybind11::hasattr(pySpatialParameters_, "extrusionFactorAtPos"))
return pySpatialParameters_.attr("extrusionFactorAtPos")(scv.dofPosition()).template cast<Scalar>();
}
// default
return 1.0;
}
/*!
* \brief Return the temperature in the given sub-control volume.
*/
template<class ElementSolution>
Scalar temperature(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
if (pySpatialParameters_)
{
if (pybind11::hasattr(pySpatialParameters_, "temperature"))
return pySpatialParameters_.attr("temperature")(element, scv, elemSol).template cast<Scalar>();
else if (pybind11::hasattr(pySpatialParameters_, "temperatureAtPos"))
return pySpatialParameters_.attr("temperatureAtPos")(scv.dofPosition()).template cast<Scalar>();
}
// default
return 283.15;
}
/*!
* \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$.
*
* The default behaviour is a constant gravity vector;
* if the <tt>Problem.EnableGravity</tt> parameter is true,
* \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$,
* else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$.
*
* \param pos the spatial position at which to evaulate the gravity vector
*/
const GravityVector& gravity(const GlobalPosition& pos) const
{ return gravity_; }
//! The finite volume grid geometry
const GridGeometry& gridGeometry() const
{ return *gridGeometry_; }
private:
std::shared_ptr<const GridGeometry> gridGeometry_;
pybind11::object pySpatialParameters_;
GravityVector gravity_;
};
// Python wrapper for the above FVProblem C++ class
template<class SpatialParams, class... options>
void registerFVSpatialParams(pybind11::handle scope, pybind11::class_<SpatialParams, options...> cls)
{
using pybind11::operator""_a;
using GridGeometry = typename SpatialParams::GridGeometry;
cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object p){
return std::make_shared<SpatialParams>(gridGeometry, p);
}));
cls.def("extrusionFactor", &SpatialParams::template extrusionFactor<decltype(std::ignore)>);
cls.def("temperature", &SpatialParams::template temperature<decltype(std::ignore)>);
cls.def("gravity", &SpatialParams::gravity);
cls.def_property_readonly("gridGeometry", &SpatialParams::gridGeometry);
}
} // namespace Dumux::Python
#endif
add_subdirectory(components)
add_subdirectory(fluidsystems)
add_subdirectory(spatialparams)
file(GLOB DUMUX_PYTHON_MATERIAL_SPATIALPARAMS_HEADERS *.hh *.inc)
install(FILES ${DUMUX_PYTHON_MATERIAL_SPATIALPARAMS_HEADERS}
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/material/spatialparams)
......@@ -35,11 +35,11 @@ namespace Dumux::Python {
* \ingroup PorousmediumflowModels
* \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>
template<class GridGeometry_, class SpatialParams_, class PrimaryVariables, bool enableInternalDirichletConstraints>
class PorousMediumFlowProblem
: public Dumux::Python::FVProblem<GridGeometry_, SpatialParams_, PrimaryVariables, enableInternalDirichletConstraints>
{
using ParentType = FVProblem<GridGeometry_, PrimaryVariables, enableInternalDirichletConstraints>;
using ParentType = Dumux::Python::FVProblem<GridGeometry_, SpatialParams_, PrimaryVariables, enableInternalDirichletConstraints>;
public:
using GridGeometry = GridGeometry_;
using SpatialParams = SpatialParams_;
......@@ -58,7 +58,7 @@ public:
PorousMediumFlowProblem(std::shared_ptr<const GridGeometry> gridGeometry,
std::shared_ptr<const SpatialParams> spatialParams,
pybind11::object pyProblem)
: ParentType(gridGeometry, pyProblem)
: ParentType(gridGeometry, spatialParams, pyProblem)
, spatialParams_(spatialParams)
{}
......@@ -79,7 +79,7 @@ void registerPorousMediumFlowProblem(pybind11::handle scope, pybind11::class_<Pr
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,
std::shared_ptr<SpatialParams> spatialParams,
pybind11::object p){
return std::make_shared<Problem>(gridGeometry, spatialParams, p);
}));
......@@ -109,7 +109,6 @@ void registerPorousMediumFlowProblem(pybind11::handle scope, pybind11::class_<Pr
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);
}
......
......@@ -21,33 +21,32 @@
* \brief TODO: docme!
*/
#ifndef DUMUX_PYTHON_MATERIAL_SPATIAL_PARAMS_HH
#define DUMUX_PYTHON_MATERIAL_SPATIAL_PARAMS_HH
#ifndef DUMUX_PYTHON_POROUSMEDIUMFLOW_FVSPATIALPARAMS_1P_HH
#define DUMUX_PYTHON_POROUSMEDIUMFLOW_FVSPATIALPARAMS_1P_HH
#include <dune/python/pybind11/pybind11.h>
#include <dune/python/pybind11/stl.h>
#include <dumux/material/spatialparams/fv1p.hh>
#include <dumux/python/common/fvspatialparams.hh>
namespace Dumux::Python {
template<class GridGeometry, class Scalar, class PT>
template<class GridGeometry_, class PT>
class FVSpatialParamsOneP
: public Dumux::FVSpatialParamsOneP<GridGeometry, Scalar, FVSpatialParamsOneP<GridGeometry, Scalar, PT>>
: public Dumux::Python::FVSpatialParams<GridGeometry_>
{
using ThisType = FVSpatialParamsOneP<GridGeometry, Scalar, PT>;
using ParentType = Dumux::FVSpatialParamsOneP<GridGeometry, Scalar, ThisType>;
using ThisType = Dumux::Python::FVSpatialParamsOneP<GridGeometry_, PT>;
using ParentType = Dumux::Python::FVSpatialParams<GridGeometry_>;
public:
using GridGeometry = GridGeometry_;
using Scalar = typename GridGeometry::GridView::ctype;
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)
: ParentType(gridGeometry, pySpatialParams)
, pySpatialParams_(pySpatialParams)
{}
......@@ -64,8 +63,8 @@ public:
template<class ElementSolution>
Scalar porosity(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol) const
const SubControlVolume& scv,
const ElementSolution& elemSol) const
{
if (pybind11::hasattr(pySpatialParams_, "porosity"))
return pySpatialParams_.attr("porosity")(element, scv, elemSol).template cast<Scalar>();
......@@ -73,24 +72,41 @@ public:
return pySpatialParams_.attr("porosityAtPos")(scv.center()).template cast<Scalar>();
}
template<class SolidSystem, class ElementSolution>
Scalar inertVolumeFraction(const Element& element,
const SubControlVolume& scv,
const ElementSolution& elemSol,
int compIdx) const
{
if (pybind11::hasattr(pySpatialParams_, "inertVolumeFraction"))
return pySpatialParams_.attr("inertVolumeFraction")(element, scv, elemSol, compIdx).template cast<Scalar>();
else if (pybind11::hasattr(pySpatialParams_, "inertVolumeFractionAtPos"))
return pySpatialParams_.attr("inertVolumeFractionAtPos")(scv.center(), compIdx).template cast<Scalar>();
else
return 1.0 - this->porosity(element, scv, elemSol);
}
static constexpr bool evaluatePermeabilityAtScvfIP()
{ return false; }
private:
pybind11::object pySpatialParams_;
};
template <class SP, class... options>
void registerOnePSpatialParams(pybind11::handle scope,
pybind11::class_<SP, options...> cls)
template <class SpatialParams, class... options>
void registerFVSpatialParamsOneP(pybind11::handle scope, pybind11::class_<SpatialParams, options...> cls)
{
using pybind11::operator""_a;
using GridGeometry = std::decay_t<decltype(std::declval<SP>().gridGeometry())>;
using GridGeometry = typename SpatialParams::GridGeometry;
cls.def(pybind11::init([](std::shared_ptr<GridGeometry> gridGeometry, pybind11::object sp){
return std::make_shared<SP>(gridGeometry, sp);
cls.def(pybind11::init([](std::shared_ptr<const GridGeometry> gridGeometry, pybind11::object p){
return std::make_shared<SpatialParams>(gridGeometry, p);
}));
cls.def("gravity", &SP::gravity);
cls.def("permeability", &SpatialParams::template permeability<decltype(std::ignore)>);
cls.def("porosity", &SpatialParams::template porosity<decltype(std::ignore)>);
cls.def("inertVolumeFraction", &SpatialParams::template inertVolumeFraction<decltype(std::ignore), decltype(std::ignore)>);
}
} // namespace Dumux::Python
......
add_python_targets(common
__init__
boundarytypes
fvproblem
fvspatialparams
properties
)
......
......@@ -7,67 +7,9 @@ from dune.generator.generator import SimpleGenerator
from dune.common.hashit import hashIt
from dumux.common.properties import Model, Property
from dumux.common.boundarytypes import BoundaryTypes
from dumux.common.fvproblem import FVProblem
from dumux.common.fvspatialparams import FVSpatialParams
from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
from ._common import *
@cppWrapperCreator
def _createFVProblemDecorator(gridGeometry, enableInternalDirichletConstraints=False):
"""A problem decorator generator for Python problems
from dumux.common import FVProblem
@FVProblem(gridGeometry)
class MyProblem:
...
"""
def createModule(numEq):
priVarType = f"Dune::FieldVector<double, {numEq}>"
ggType = gridGeometry._typeName
enableIntDirConstraint = "true" if enableInternalDirichletConstraints else "false"
problemType = f"Dumux::Python::FVProblem<{ggType}, {priVarType}, {enableIntDirConstraint}>"
includes = gridGeometry._includes + ["dumux/python/common/fvproblem.hh"]
moduleName = "fvproblem_" + hashIt(problemType)
holderType = f"std::shared_ptr<{problemType}>"
generator = SimpleGenerator("FVProblem", "Dumux::Python")
module = generator.load(includes, problemType, moduleName, options=[holderType])
return module
def decorateFVProblem(cls):
module = createModule(cls.numEq)
def createFVProblem():
return module.FVProblem(gridGeometry, cls())
return createFVProblem
return decorateFVProblem
@cppWrapperClassAlias(creator=_createFVProblemDecorator)
class FVProblem:
"""Class alias used to decorate a Python finite volume problem"""
@cppWrapperCreator
def _createBoundaryTypes(numEq=1):
"""Create BoundaryTypes instances"""
# only compile this once per numEq
cacheKey = f"BoundaryTypes_{numEq}"
try:
return globals()[cacheKey]()
except KeyError:
includes = ["dumux/python/common/boundarytypes.hh"]
typeName = f"Dumux::BoundaryTypes<{numEq}>"
moduleName = "boundarytypes_" + hashIt(typeName)
generator = SimpleGenerator("BoundaryTypes", "Dumux::Python")
module = generator.load(includes, typeName, moduleName)
globals().update({cacheKey: module.BoundaryTypes})
return globals()[cacheKey]()
@cppWrapperClassAlias(creator=_createBoundaryTypes)
class BoundaryTypes:
"""Class alias used to create a BoundaryTypes instance"""
"""
Boundary types generator
"""
from dune.generator.generator import SimpleGenerator
from dune.common.hashit import hashIt
from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
@cppWrapperCreator
def _createBoundaryTypes(numEq=1):
"""Create BoundaryTypes instances"""
# only compile this once per numEq
cacheKey = f"BoundaryTypes_{numEq}"
try:
return globals()[cacheKey]()
except KeyError:
includes = ["dumux/python/common/boundarytypes.hh"]
typeName = f"Dumux::BoundaryTypes<{numEq}>"
moduleName = "boundarytypes_" + hashIt(typeName)
generator = SimpleGenerator("BoundaryTypes", "Dumux::Python")
module = generator.load(includes, typeName, moduleName)
globals().update({cacheKey: module.BoundaryTypes})
return globals()[cacheKey]()
@cppWrapperClassAlias(creator=_createBoundaryTypes)
class BoundaryTypes:
"""Class alias used to create a BoundaryTypes instance"""
"""
Finite volume problem generator
"""
from dune.generator.generator import SimpleGenerator
from dune.common.hashit import hashIt
from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias
@cppWrapperCreator
def _createFVProblemDecorator(
gridGeometry, spatialParams=None, enableInternalDirichletConstraints=False
):
"""A problem decorator generator for Python problems
from dumux.common import FVProblem
@FVProblem(gridGeometry)
class MyProblem:
...
"""
def createModule(numEq):
ggType = gridGeometry._typeName
if spatialParams is not None:
spType = spatialParams._typeName