From f97a945f34b0a6aef051c05b864bd5f666ad7097 Mon Sep 17 00:00:00 2001 From: Timo Koch <timokoch@math.uio.no> Date: Mon, 2 May 2022 15:14:20 +0000 Subject: [PATCH] Merge branch 'feature/spatialparams_python' into 'master' Spatial parameters for Python bindings Closes #1145 See merge request dumux-repositories/dumux!3069 (cherry picked from commit 790bdaeeea272a72c94facb388db0f462a54b498) 8e464d6c [dumux][python][spatialparams] update the problem and spatialparams classes 2e25f728 [python][material] remove old spatial params files 56a3c69b [python] Make spatialparams optional and make separate modules in python/common 97d568bb [python] Make spatial params const again --- dumux/python/common/fvproblem.hh | 58 ++++-- dumux/python/common/fvspatialparams.hh | 166 ++++++++++++++++++ dumux/python/material/CMakeLists.txt | 1 - .../material/spatialparams/CMakeLists.txt | 3 - dumux/python/porousmediumflow/problem.hh | 13 +- .../spatialparams.hh | 62 ++++--- python/dumux/common/CMakeLists.txt | 3 + python/dumux/common/__init__.py | 64 +------ python/dumux/common/boundarytypes.py | 31 ++++ python/dumux/common/fvproblem.py | 58 ++++++ python/dumux/common/fvspatialparams.py | 44 +++++ python/dumux/material/CMakeLists.txt | 1 - python/dumux/material/__init__.py | 1 - .../material/spatialparams/CMakeLists.txt | 3 - .../dumux/material/spatialparams/__init__.py | 77 -------- python/dumux/porousmediumflow/__init__.py | 74 +++++++- test/python/test_1p.py | 20 +-- test/python/test_explicit_transport_cctpfa.py | 8 +- test/python/test_fvproblem.py | 2 +- 19 files changed, 478 insertions(+), 211 deletions(-) create mode 100644 dumux/python/common/fvspatialparams.hh delete mode 100644 dumux/python/material/spatialparams/CMakeLists.txt rename dumux/python/{material/spatialparams => porousmediumflow}/spatialparams.hh (57%) create mode 100644 python/dumux/common/boundarytypes.py create mode 100644 python/dumux/common/fvproblem.py create mode 100644 python/dumux/common/fvspatialparams.py delete mode 100644 python/dumux/material/spatialparams/CMakeLists.txt delete mode 100644 python/dumux/material/spatialparams/__init__.py diff --git a/dumux/python/common/fvproblem.hh b/dumux/python/common/fvproblem.hh index 86eb699bc1..5da2a9f348 100644 --- a/dumux/python/common/fvproblem.hh +++ b/dumux/python/common/fvproblem.hh @@ -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 diff --git a/dumux/python/common/fvspatialparams.hh b/dumux/python/common/fvspatialparams.hh new file mode 100644 index 0000000000..9b69fe339b --- /dev/null +++ b/dumux/python/common/fvspatialparams.hh @@ -0,0 +1,166 @@ +// -*- 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 diff --git a/dumux/python/material/CMakeLists.txt b/dumux/python/material/CMakeLists.txt index f1c168afc3..8400d5f3eb 100644 --- a/dumux/python/material/CMakeLists.txt +++ b/dumux/python/material/CMakeLists.txt @@ -1,3 +1,2 @@ add_subdirectory(components) add_subdirectory(fluidsystems) -add_subdirectory(spatialparams) diff --git a/dumux/python/material/spatialparams/CMakeLists.txt b/dumux/python/material/spatialparams/CMakeLists.txt deleted file mode 100644 index 08816c5b1a..0000000000 --- a/dumux/python/material/spatialparams/CMakeLists.txt +++ /dev/null @@ -1,3 +0,0 @@ -file(GLOB DUMUX_PYTHON_MATERIAL_SPATIALPARAMS_HEADERS *.hh *.inc) -install(FILES ${DUMUX_PYTHON_MATERIAL_SPATIALPARAMS_HEADERS} - DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/python/material/spatialparams) diff --git a/dumux/python/porousmediumflow/problem.hh b/dumux/python/porousmediumflow/problem.hh index 7b369fade4..d97b8a28eb 100644 --- a/dumux/python/porousmediumflow/problem.hh +++ b/dumux/python/porousmediumflow/problem.hh @@ -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); } diff --git a/dumux/python/material/spatialparams/spatialparams.hh b/dumux/python/porousmediumflow/spatialparams.hh similarity index 57% rename from dumux/python/material/spatialparams/spatialparams.hh rename to dumux/python/porousmediumflow/spatialparams.hh index 448a981491..2199309279 100644 --- a/dumux/python/material/spatialparams/spatialparams.hh +++ b/dumux/python/porousmediumflow/spatialparams.hh @@ -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 diff --git a/python/dumux/common/CMakeLists.txt b/python/dumux/common/CMakeLists.txt index 21bad44ccf..31a5167f3d 100644 --- a/python/dumux/common/CMakeLists.txt +++ b/python/dumux/common/CMakeLists.txt @@ -1,5 +1,8 @@ add_python_targets(common __init__ + boundarytypes + fvproblem + fvspatialparams properties ) diff --git a/python/dumux/common/__init__.py b/python/dumux/common/__init__.py index 879a0f7a8f..f2fa10c161 100644 --- a/python/dumux/common/__init__.py +++ b/python/dumux/common/__init__.py @@ -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""" diff --git a/python/dumux/common/boundarytypes.py b/python/dumux/common/boundarytypes.py new file mode 100644 index 0000000000..7029b5ca2c --- /dev/null +++ b/python/dumux/common/boundarytypes.py @@ -0,0 +1,31 @@ +""" +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""" diff --git a/python/dumux/common/fvproblem.py b/python/dumux/common/fvproblem.py new file mode 100644 index 0000000000..416500c316 --- /dev/null +++ b/python/dumux/common/fvproblem.py @@ -0,0 +1,58 @@ +""" +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 + else: + spType = f"Dumux::Python::FVSpatialParams<{ggType}>" + + priVarType = f"Dune::FieldVector<double, {numEq}>" + enableIntDirConstraint = "true" if enableInternalDirichletConstraints else "false" + problemType = ( + f"Dumux::Python::FVProblem<{ggType}, {spType}, {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(): + if spatialParams is not None: + return module.FVProblem(gridGeometry, spatialParams, cls()) + + return module.FVProblem(gridGeometry, cls()) + + return createFVProblem + + return decorateFVProblem + + +@cppWrapperClassAlias(creator=_createFVProblemDecorator) +class FVProblem: + """Class alias used to decorate a Python finite volume problem""" diff --git a/python/dumux/common/fvspatialparams.py b/python/dumux/common/fvspatialparams.py new file mode 100644 index 0000000000..812df3f3b6 --- /dev/null +++ b/python/dumux/common/fvspatialparams.py @@ -0,0 +1,44 @@ +""" +Finite volume spatial params generator +""" + +from dune.generator.generator import SimpleGenerator +from dune.common.hashit import hashIt + +from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias + + +@cppWrapperCreator +def _createFVSpatialParamsDecorator(gridGeometry): + """A spatial params decorator generator for Python spatial params + + from dumux.common import FVSpatialParams + @FVSpatialParams(gridGeometry) + class MySpatialParams: + ... + """ + + def createModule(): + ggType = gridGeometry._typeName + spatialParamsType = f"Dumux::Python::FVSpatialParams<{ggType}>" + includes = gridGeometry._includes + ["dumux/python/common/fvspatialparams.hh"] + moduleName = "fvspatialparams_" + hashIt(spatialParamsType) + holderType = f"std::shared_ptr<{spatialParamsType}>" + generator = SimpleGenerator("FVSpatialParams", "Dumux::Python") + module = generator.load(includes, spatialParamsType, moduleName, options=[holderType]) + return module + + def decorateFVSpatialParams(cls): + module = createModule() + + def createFVSpatialParams(): + return module.FVSpatialParams(gridGeometry, cls()) + + return createFVSpatialParams + + return decorateFVSpatialParams + + +@cppWrapperClassAlias(creator=_createFVSpatialParamsDecorator) +class FVSpatialParams: + """Class alias used to decorate a Python Finite Volume Spatial Params""" diff --git a/python/dumux/material/CMakeLists.txt b/python/dumux/material/CMakeLists.txt index c6683fdb0b..71b18100e7 100644 --- a/python/dumux/material/CMakeLists.txt +++ b/python/dumux/material/CMakeLists.txt @@ -1,6 +1,5 @@ add_subdirectory(components) add_subdirectory(fluidsystems) -add_subdirectory(spatialparams) add_python_targets(material __init__ diff --git a/python/dumux/material/__init__.py b/python/dumux/material/__init__.py index 01534f2422..f3a972ef3f 100644 --- a/python/dumux/material/__init__.py +++ b/python/dumux/material/__init__.py @@ -2,4 +2,3 @@ from dumux.material.fluidsystems import FluidSystem from dumux.material.components import Component -from dumux.material.spatialparams import OnePSpatialParams diff --git a/python/dumux/material/spatialparams/CMakeLists.txt b/python/dumux/material/spatialparams/CMakeLists.txt deleted file mode 100644 index db81731a16..0000000000 --- a/python/dumux/material/spatialparams/CMakeLists.txt +++ /dev/null @@ -1,3 +0,0 @@ -add_python_targets(spatialparams - __init__ -) diff --git a/python/dumux/material/spatialparams/__init__.py b/python/dumux/material/spatialparams/__init__.py deleted file mode 100644 index 49c9504581..0000000000 --- a/python/dumux/material/spatialparams/__init__.py +++ /dev/null @@ -1,77 +0,0 @@ -"""Spatial parameter classes""" - -import numpy as np -from dune.generator.generator import SimpleGenerator -from dune.common.hashit import hashIt -from dumux.common import Property -from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias - - -@cppWrapperCreator -def _createOnePSpatialParamsDecorator(*, gridGeometry, scalar: Property = None): - """Turn a Python spatial parameter class into an C++/Python hybrid class""" - - dim = gridGeometry.gridView.dimensionworld - includes = gridGeometry._includes + [ - "dumux/python/material/spatialparams/spatialparams.hh", - "dune/common/fmatrix.hh", - ] - - if scalar is None: - scalar = Property.fromCppType("double") - scalarType = scalar.cppType - - permType = f"Dune::FieldMatrix<{scalarType}, {dim}, {dim}>" - typeName = ( - f"Dumux::Python::FVSpatialParamsOneP<{gridGeometry._typeName}, {scalarType}, {permType}>" - ) - moduleName = f"spatialparams_{hashIt(typeName)}" - - def decorateOnePSpatialParams(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() - - return permeabilityValue - - class Permeability: - """Permeability decorator to make sure permeability has correct type""" - - def __init__(self, permeabilityFunction): - self.permeabilityFunction = permeabilityFunction - - def __call__(self, element, scv, elemSol): - result = self.permeabilityFunction(element, scv, elemSol) - return maybeConvertScalarToMatrix(result) - - class PermeabilityAtPos: - """PermeabilityAtPos decorator to make sure permeability has correct type""" - - def __init__(self, permeabilityFunction): - self.permeabilityFunction = permeabilityFunction - - def __call__(self, globalPos): - result = self.permeabilityFunction(globalPos) - return maybeConvertScalarToMatrix(result) - - def createSpatialParams(): - spatialParams = cls() - if hasattr(cls, "permeability"): - cls.permeability = Permeability(spatialParams.permeability) - if hasattr(cls, "permeabilityAtPos"): - cls.permeabilityAtPos = PermeabilityAtPos(spatialParams.permeabilityAtPos) - return module.OnePSpatialParams(gridGeometry, spatialParams) - - return createSpatialParams - - return decorateOnePSpatialParams - - -@cppWrapperClassAlias(creator=_createOnePSpatialParamsDecorator) -class OnePSpatialParams: - """Class alias used to decorate Python spatial parameter implementations""" diff --git a/python/dumux/porousmediumflow/__init__.py b/python/dumux/porousmediumflow/__init__.py index 3a4c9c6729..42892d308b 100644 --- a/python/dumux/porousmediumflow/__init__.py +++ b/python/dumux/porousmediumflow/__init__.py @@ -1,7 +1,9 @@ """Classes and functions related to the porousmedium flow models""" +import numpy as np from dune.generator.generator import SimpleGenerator from dune.common.hashit import hashIt +from dumux.common import Property from dumux.wrapping import cppWrapperCreator, cppWrapperClassAlias @@ -24,7 +26,7 @@ def _createPorousMediumFlowProblemDecorator( spType = spatialParams._typeName enableIDC = "true" if enableInternalDirichletConstraints else "false" problemType = ( - "Dumux::Python::PorousMediumFlowProblem" f"<{ggType}, {priVars}, {spType}, {enableIDC}>" + "Dumux::Python::PorousMediumFlowProblem" f"<{ggType}, {spType}, {priVars}, {enableIDC}>" ) includes = [ *(gridGeometry._includes), @@ -79,3 +81,73 @@ def _createPorousMediumFlowVelocityOutput(*, gridVariables): @cppWrapperClassAlias(creator=_createPorousMediumFlowVelocityOutput) class PorousMediumFlowVelocityOutput: """A class alias used to create PorousMediumFlowVelocityOutput instances""" + + +@cppWrapperCreator +def _createFVSpatialParamsOnePDecorator(gridGeometry): + """A spatial params decorator generator for Python spatial params + + from dumux.common import FVSpatialParamsOneP + FVSpatialParamsOneP(gridGeometry) + class MySpatialParams: + ... + """ + dim = gridGeometry.gridView.dimensionworld + scalar = Property.fromCppType("double") + scalarType = scalar.cppType + permType = f"Dune::FieldMatrix<{scalarType}, {dim}, {dim}>" + + def decorateFVSpatialParamsOneP(cls): + ggType = gridGeometry._typeName + spatialParamsType = f"Dumux::Python::FVSpatialParamsOneP<{ggType}, {permType}>" + includes = gridGeometry._includes + ["dumux/python/porousmediumflow/spatialparams.hh"] + includes += ["dumux/python/common/fvspatialparams.hh"] + moduleName = "fvspatialparamsonep_" + hashIt(spatialParamsType) + holderType = f"std::shared_ptr<{spatialParamsType}>" + generator = SimpleGenerator("FVSpatialParamsOneP", "Dumux::Python") + module = generator.load(includes, spatialParamsType, moduleName, options=[holderType]) + + def maybeConvertScalarToMatrix(permeabilityValue): + if isinstance(permeabilityValue, float): + matrix = np.zeros(shape=(dim, dim)) + np.fill_diagonal(matrix, permeabilityValue) + return matrix.tolist() + + return permeabilityValue + + class Permeability: + """Permeability decorator to make sure permeability has correct type""" + + def __init__(self, permeabilityFunction): + self.permeabilityFunction = permeabilityFunction + + def __call__(self, element, scv, elemSol): + result = self.permeabilityFunction(element, scv, elemSol) + return maybeConvertScalarToMatrix(result) + + class PermeabilityAtPos: + """PermeabilityAtPos decorator to make sure permeability has correct type""" + + def __init__(self, permeabilityFunction): + self.permeabilityFunction = permeabilityFunction + + def __call__(self, globalPos): + result = self.permeabilityFunction(globalPos) + return maybeConvertScalarToMatrix(result) + + def createFVSpatialParamsOneP(): + spatialParams = cls() + if hasattr(cls, "permeability"): + cls.permeability = Permeability(spatialParams.permeability) + if hasattr(cls, "permeabilityAtPos"): + cls.permeabilityAtPos = PermeabilityAtPos(spatialParams.permeabilityAtPos) + return module.FVSpatialParamsOneP(gridGeometry, cls()) + + return createFVSpatialParamsOneP + + return decorateFVSpatialParamsOneP + + +@cppWrapperClassAlias(creator=_createFVSpatialParamsOnePDecorator) +class FVSpatialParamsOneP: + """Class alias used to decorate a Python Finite Volume Spatial Params""" diff --git a/test/python/test_1p.py b/test/python/test_1p.py index 89b14accbc..de0c1d0037 100755 --- a/test/python/test_1p.py +++ b/test/python/test_1p.py @@ -9,8 +9,12 @@ from dumux.common import initParameters, printParameters, getParam from dumux.common import BoundaryTypes, Model, Property from dumux.discretization import GridGeometry, GridVariables from dumux.assembly import FVAssembler -from dumux.porousmediumflow import PorousMediumFlowProblem, PorousMediumFlowVelocityOutput -from dumux.material import FluidSystem, Component, OnePSpatialParams +from dumux.porousmediumflow import ( + PorousMediumFlowProblem, + PorousMediumFlowVelocityOutput, + FVSpatialParamsOneP, +) +from dumux.material import FluidSystem, Component from dumux.io import VtkOutputModule # Initialize parameters @@ -52,7 +56,7 @@ model["FluidSystem"] = Property.fromInstance(onePLiquid) # Define the spatial parameters -@OnePSpatialParams(gridGeometry=gridGeometry) +@FVSpatialParamsOneP(gridGeometry=gridGeometry) class SpatialParams: dimWorld = gridGeometry.gridView.dimWorld lensLowerLeft = [0.2, 0.2] @@ -80,13 +84,11 @@ class SpatialParams: return 0.4 -# and set them as a model property spatialParams = SpatialParams() model["SpatialParams"] = Property.fromInstance(spatialParams) - # Define the problem -@PorousMediumFlowProblem(gridGeometry, spatialParams) +@PorousMediumFlowProblem(gridGeometry=gridGeometry, spatialParams=spatialParams) class Problem: numEq = 1 @@ -108,12 +110,6 @@ class Problem: def sourceAtPos(self, globalPos): return 0.0 - def extrusionFactor(self, element, scv): - return 1.0 - - def temperatureAtPos(self, globalPos): - return 283.15 - def name(self): return "python_problem" diff --git a/test/python/test_explicit_transport_cctpfa.py b/test/python/test_explicit_transport_cctpfa.py index 4ad605c3a9..670716406a 100755 --- a/test/python/test_explicit_transport_cctpfa.py +++ b/test/python/test_explicit_transport_cctpfa.py @@ -29,12 +29,12 @@ gridGeometry = GridGeometry(gridView, discMethod="cctpfa") elementMapper = gridView.indexSet -############################################## -# Define problem (inital/boundary condtions) # -############################################## +################################################ +# Define problem (initial/boundary conditions) # +################################################ -@FVProblem(gridGeometry) +@FVProblem(gridGeometry=gridGeometry) class Problem: numEq = 1 diff --git a/test/python/test_fvproblem.py b/test/python/test_fvproblem.py index 12adceaad7..526a23aab2 100755 --- a/test/python/test_fvproblem.py +++ b/test/python/test_fvproblem.py @@ -27,7 +27,7 @@ gridView = structuredGrid([0, 0, 0], [1, 1, 1], [3, 3, 3]) gridGeometry = GridGeometry(gridView, discMethod="box") -@FVProblem(gridGeometry) +@FVProblem(gridGeometry=gridGeometry) class Problem: numEq = 2 -- GitLab