diff --git a/dumux/python/common/fvproblem.hh b/dumux/python/common/fvproblem.hh index c654b34f0f4c8163eb9a6f92b51e06ee51d12005..bf8aaba7bf0d6af994d1c539594a730d3afed192 100644 --- a/dumux/python/common/fvproblem.hh +++ b/dumux/python/common/fvproblem.hh @@ -28,9 +28,6 @@ #include <memory> #include <tuple> -#include <iostream> // for debug stuff below -#include <dune/python/pybind11/iostream.h> // for debug stuff below - #include <dune/common/fvector.hh> #include <dune/common/exceptions.hh> #include <dune/python/pybind11/pybind11.h> @@ -193,64 +190,6 @@ void registerFVProblem(pybind11::handle scope, pybind11::class_<Problem, options cls.def("gridGeometry", &Problem::gridGeometry); } - -///////////////////////////////////////////////////////// -/// Some debugging/testing stuff -/////////////////////////////////////////////////////////// -template<class Problem> -void printBoundaryStuff(const Problem& problem) -{ - const auto& gg = problem.gridGeometry(); - for (const auto& element : elements(gg.gridView())) - { - auto fvGeometry = localView(gg); - fvGeometry.bindElement(element); - - for (const auto& scv : scvs(fvGeometry)) - { - const auto boundaryTypes = problem.boundaryTypes(element, scv); - std::cout << "-- scv at " << scv.dofPosition() << ": " - << " isNeumann: " << std::boolalpha << boundaryTypes.hasNeumann() - << " isDirichlet: " << boundaryTypes.hasDirichlet() << std::endl; - - if (boundaryTypes.hasDirichlet()) - std::cout << " -- Dirichlet values: " << problem.dirichlet(element, scv) << std::endl; - } - } -} - -template<class P> -class PrintBoundaryStuff -{ -public: - using Problem = P; - - PrintBoundaryStuff(std::shared_ptr<const Problem> problem) - : problem_(problem) - {} - - void print() - { - printBoundaryStuff(*problem_); - } - -private: - std::shared_ptr<const Problem> problem_; -}; - -template<class T, class... options> -void registerPrintBoundaryStuff(pybind11::handle scope, pybind11::class_<T, options...> cls) -{ - using Problem = typename T::Problem; - cls.def(pybind11::init([](std::shared_ptr<const Problem> problem){ - return std::make_unique<T>(problem); - })); - cls.def("print", [](T& self){ - pybind11::scoped_ostream_redirect stream(std::cout, pybind11::module::import("sys").attr("stdout")); - self.print(); - }); -} - } // end namespace Dumux::Python #endif diff --git a/python/dumux/common/__init__.py b/python/dumux/common/__init__.py index 4d9429121c48cc88f85561e22da3bf1ce54f7bd6..23cdbfa5710456ff954c9281e27e7df725cb6ff4 100644 --- a/python/dumux/common/__init__.py +++ b/python/dumux/common/__init__.py @@ -46,13 +46,3 @@ def BoundaryTypes(numEq=1): module = generator.load(includes, typeName, moduleName) globals().update({cacheKey : module.BoundaryTypes}) return globals()[cacheKey]() - - -# debugging/testing code -def PrintBoundaryStuff(problem): - includes = problem._includes + ["dumux/python/common/fvproblem.hh"] - typeName = "Dumux::Python::PrintBoundaryStuff<{}>".format(problem._typeName) - moduleName = moduleName = "printbs_" + hashIt(problem._typeName) - generator = SimpleGenerator("PrintBoundaryStuff", "Dumux::Python") - module = generator.load(includes, typeName, moduleName) - return module.PrintBoundaryStuff(problem) diff --git a/test/python/test_boundaryloop.hh b/test/python/test_boundaryloop.hh new file mode 100644 index 0000000000000000000000000000000000000000..9642f390b8843bf3859381009e8eb53aa130d40a --- /dev/null +++ b/test/python/test_boundaryloop.hh @@ -0,0 +1,110 @@ +// -*- 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_TEST_PYTHON_FVPROBLEM_BOUNDARY_LOOP_HH +#define DUMUX_TEST_PYTHON_FVPROBLEM_BOUNDARY_LOOP_HH + +#include <string> +#include <memory> +#include <tuple> + +#include <iostream> +#include <iomanip> +#include <dune/python/pybind11/iostream.h> + +#include <dune/common/fvector.hh> +#include <dune/common/exceptions.hh> +#include <dune/python/pybind11/pybind11.h> + +#include <dumux/common/boundarytypes.hh> +#include <dumux/python/common/boundarytypes.hh> + +namespace Dumux::Python { + +///////////////////////////////////////////////////////// +/// Some debugging/testing stuff +/////////////////////////////////////////////////////////// +template<class Problem> +void printProblemTest(const Problem& problem) +{ + std::size_t numNeumann = 0; + std::size_t numDirichlet = 0; + double totalSource = 0; + const auto& gg = problem.gridGeometry(); + for (const auto& element : elements(gg.gridView())) + { + auto fvGeometry = localView(gg); + fvGeometry.bindElement(element); + + for (const auto& scv : scvs(fvGeometry)) + { + const auto boundaryTypes = problem.boundaryTypes(element, scv); + + if (boundaryTypes.hasDirichlet()) + numDirichlet += 1; + else if (boundaryTypes.hasNeumann()) + numNeumann += 1; + + totalSource += problem.sourceAtPos(scv.center())[0]*scv.volume(); + } + } + + std::cout << "[cpp] Found " << numNeumann << " Neumann faces and " << numDirichlet << " Dirichlet faces" << std::endl; + std::cout << "[cpp] Total source " << totalSource << std::fixed << std::setprecision(3) << " kg/s" << std::endl; +} + +template<class P> +class PrintProblemTest +{ +public: + using Problem = P; + + PrintProblemTest(std::shared_ptr<const Problem> problem) + : problem_(problem) + {} + + void print() + { + printProblemTest(*problem_); + } + +private: + std::shared_ptr<const Problem> problem_; +}; + +template<class T, class... options> +void registerPrintProblemTest(pybind11::handle scope, pybind11::class_<T, options...> cls) +{ + using Problem = typename T::Problem; + cls.def(pybind11::init([](std::shared_ptr<const Problem> problem){ + return std::make_unique<T>(problem); + })); + cls.def("print", [](T& self){ + pybind11::scoped_ostream_redirect stream(std::cout, pybind11::module::import("sys").attr("stdout")); + self.print(); + }); +} + +} // end namespace Dumux::Python + +#endif diff --git a/test/python/test_py_fvproblem.py b/test/python/test_py_fvproblem.py index f676692d0b10b94dc7ef9e79e6139c8df3e1ba1d..74ecbf10570cf3fd7fae73c9178b5abc47e5d448 100644 --- a/test/python/test_py_fvproblem.py +++ b/test/python/test_py_fvproblem.py @@ -1,9 +1,25 @@ #!/usr/bin/env python3 -import time +# Compile some C++ code for reference +from dumux.common import * +from dune.generator.generator import SimpleGenerator +from dune.common.hashit import hashIt + +# debugging/testing code +def PrintProblemTest(problem): + includes = problem._includes + ["test/python/test_boundaryloop.hh"] + typeName = "Dumux::Python::PrintProblemTest<{}>".format(problem._typeName) + moduleName = moduleName = "printbs_" + hashIt(problem._typeName) + generator = SimpleGenerator("PrintProblemTest", "Dumux::Python") + module = generator.load(includes, typeName, moduleName) + return module.PrintProblemTest(problem) + +############################################################ +# The actual Python test +############################################################ from dune.grid import structuredGrid from dumux.discretization import GridGeometry -from dumux.common import BoundaryTypes, FVProblem, PrintBoundaryStuff +from dumux.common import BoundaryTypes, FVProblem gridView = structuredGrid([0,0,0],[1,1,1],[3,3,3]) @@ -26,6 +42,9 @@ class Problem: else: return [1.0, 0.0] + def sourceAtPos(self, globalPos): + return [globalPos[0]] + problem = Problem() print("Name of the problem: {}".format(problem.name)) print("-- Number of equations: {}".format(problem.numEq)) @@ -34,16 +53,23 @@ print() # test C++/Python compatibility of the hybrid class # print some info from C++ -printer = PrintBoundaryStuff(problem) +printer = PrintProblemTest(problem) printer.print() # print the same info from Python -start = time.time() +numNeumann = 0 +numDirichlet = 0 +totalSource = 0 for e in gridView.elements: fvGeometry = problem.gridGeometry().localView() fvGeometry.bind(e) for scv in fvGeometry.scvs(): bTypes = problem.boundaryTypes(element=e, scv=scv) - print("-- scv at {}: isNeumann: {}, isDirichlet: {}".format(scv.dofPosition(), bTypes.isNeumann(), bTypes.isDirichlet())) if bTypes.isDirichlet(): - print(" -- Dirichlet values: {}".format(problem.dirichlet(element=e, scv=scv))) + numDirichlet += 1 + elif bTypes.isNeumann(): + numNeumann += 1 + totalSource += problem.sourceAtPos(scv.center())[0]*scv.volume() + +print("[python] Found {} Neumann faces and {} Dirichlet faces".format(numNeumann, numDirichlet)) +print("[python] Total source {:.2f} kg/s".format(totalSource))