From a10a415bed310d5cafa9698039c3b52b6df21df6 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 31 May 2019 23:14:46 +0200
Subject: [PATCH] [common] Add function to integrate FV-style solutions
---
dumux/common/integrate.hh | 120 +++++++++++++++++-
test/common/functions/test_function_l2norm.cc | 53 +++++---
2 files changed, 152 insertions(+), 21 deletions(-)
diff --git a/dumux/common/integrate.hh b/dumux/common/integrate.hh
index 9ea524a6da..3a68c6d83c 100644
--- a/dumux/common/integrate.hh
+++ b/dumux/common/integrate.hh
@@ -28,17 +28,117 @@
#include <type_traits>
#include <dune/geometry/quadraturerules.hh>
+#include <dune/common/concept.hh>
+
+#if HAVE_DUNE_FUNCTIONS
+#include <dune/functions/gridfunctions/gridfunction.hh>
+#endif
+
+#include <dumux/discretization/evalsolution.hh>
+#include <dumux/discretization/elementsolution.hh>
namespace Dumux {
+namespace Impl {
+
+struct HasLocalFunction
+{
+ template<class F>
+ auto require(F&& f) -> decltype(
+ localFunction(f),
+ localFunction(f).unbind()
+ );
+};
+
+template<class F>
+static constexpr bool hasLocalFunction()
+{ return Dune::models<HasLocalFunction, F>(); }
+
+} // end namespace Impl
+
+/*!
+ * \brief Integrate a grid function over a grid view
+ * \param gv the grid view
+ * \param f the grid function
+ * \param order the order of the quadrature rule
+ */
+template<class GridGeometry, class SolutionVector,
+ typename std::enable_if_t<!Impl::hasLocalFunction<SolutionVector>(), int> = 0>
+auto integrateGridFunction(const GridGeometry& gg,
+ SolutionVector&& sol,
+ std::size_t order)
+{
+ using Scalar = std::decay_t<decltype(sol[0][0])>;
+ using GridView = typename GridGeometry::GridView;
+
+ Scalar integral = 0.0;
+ for (const auto& element : elements(gg.gridView()))
+ {
+ const auto elemSol = elementSolution(element, sol, gg);
+ const auto geometry = element.geometry();
+ const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
+ for (auto&& qp : quad)
+ {
+ auto value = evalSolution(element, geometry, gg, elemSol, geometry.global(qp.position()));
+ value *= qp.weight()*geometry.integrationElement(qp.position());
+ integral += value;
+ }
+ }
+ return integral;
+}
+
+/*!
+ * \brief Integrate a function over a grid view
+ * \param gv the grid view
+ * \param f the first function
+ * \param g the second function
+ * \param order the order of the quadrature rule
+ * \note dune functions currently doesn't support composing two functions
+ */
+template<class GridGeometry, class Sol1, class Sol2,
+ typename std::enable_if_t<!Impl::hasLocalFunction<Sol1>(), int> = 0>
+auto integrateL2Error(const GridGeometry& gg,
+ const Sol1& sol1,
+ const Sol2& sol2,
+ std::size_t order)
+{
+ using Scalar = std::decay_t<decltype(sol1[0][0])>;
+ using GridView = typename GridGeometry::GridView;
+
+ Scalar l2norm = 0.0;
+ for (const auto& element : elements(gg.gridView()))
+ {
+ const auto elemSol1 = elementSolution(element, sol1, gg);
+ const auto elemSol2 = elementSolution(element, sol2, gg);
+
+ const auto geometry = element.geometry();
+ const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
+ for (auto&& qp : quad)
+ {
+ const auto& globalPos = geometry.global(qp.position());
+ const auto value1 = evalSolution(element, geometry, gg, elemSol1, globalPos);
+ const auto value2 = evalSolution(element, geometry, gg, elemSol2, globalPos);
+ const auto error = (value1 - value2);
+ l2norm += error*error*qp.weight()*geometry.integrationElement(qp.position());
+ }
+ }
+ using std::sqrt;
+ return sqrt(l2norm);
+}
+
+#if HAVE_DUNE_FUNCTIONS
/*!
* \brief Integrate a grid function over a grid view
* \param gv the grid view
* \param f the grid function
* \param order the order of the quadrature rule
+ * \note overload for a Dune::Funtions::GridFunction
*/
-template<class GridView, class F>
-auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
+template<class GridView, class F,
+ typename std::enable_if_t<Impl::hasLocalFunction<F>(), int> = 0>
+auto integrateGridFunction(const GridView& gv,
+ F&& f,
+ std::size_t order)
{
auto fLocal = localFunction(f);
@@ -54,7 +154,11 @@ auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
const auto geometry = element.geometry();
const auto& quad = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(geometry.type(), order);
for (auto&& qp : quad)
- integral += fLocal(qp.position())*qp.weight()*geometry.integrationElement(qp.position());
+ {
+ auto value = fLocal(qp.position());
+ value *= qp.weight()*geometry.integrationElement(qp.position());
+ integral += value;
+ }
fLocal.unbind();
}
@@ -67,10 +171,15 @@ auto integrateGridFunction(const GridView& gv, F&& f, std::size_t order)
* \param f the first function
* \param g the second function
* \param order the order of the quadrature rule
+ * \note overload for a Dune::Funtions::GridFunction
* \note dune functions currently doesn't support composing two functions
*/
-template<class GridView, class F, class G>
-auto integrateL2Error(const GridView& gv, F&& f, G&& g, std::size_t order)
+template<class GridView, class F, class G,
+ typename std::enable_if_t<Impl::hasLocalFunction<F>(), int> = 0>
+auto integrateL2Error(const GridView& gv,
+ F&& f,
+ G&& g,
+ std::size_t order)
{
auto fLocal = localFunction(f);
auto gLocal = localFunction(g);
@@ -102,6 +211,7 @@ auto integrateL2Error(const GridView& gv, F&& f, G&& g, std::size_t order)
using std::sqrt;
return sqrt(l2norm);
}
+#endif
} // end namespace Dumux
diff --git a/test/common/functions/test_function_l2norm.cc b/test/common/functions/test_function_l2norm.cc
index f4911dda84..baaba12f26 100644
--- a/test/common/functions/test_function_l2norm.cc
+++ b/test/common/functions/test_function_l2norm.cc
@@ -12,7 +12,8 @@
#include <dumux/common/integrate.hh>
#include <dumux/multidomain/glue.hh>
#include <dumux/discretization/projection/projector.hh>
-#include <dumux/discretization/fem/fegridgeometry.hh>
+#include <dumux/discretization/basegridgeometry.hh>
+#include <dumux/discretization/box/fvgridgeometry.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
int main (int argc, char *argv[]) try
@@ -35,38 +36,58 @@ int main (int argc, char *argv[]) try
const auto gvFine = gridManagerFine.grid().leafGridView();
const auto gvCoarse = gridManagerCoarse.grid().leafGridView();
- // a quadratic function on the grid
- auto f = [](const auto& pos){ return pos.two_norm2(); };
-
- // make discrete functions
+ // get some types
using R = Dune::FieldVector<double, 1>;
using SolutionVector = Dune::BlockVector<R>;
using GridView = Grid::LeafGridView;
using namespace Dune::Functions;
using P2 = LagrangeBasis<GridView, 2>;
- auto basisCoarse = std::make_shared<P2>(gvCoarse), basisFine = std::make_shared<P2>(gvFine);
+ // construct glue, any grid geometry should be fine here
+ Dumux::BaseGridGeometry<GridView, Dumux::DefaultMapperTraits<GridView>> ggCoarse(gvCoarse), ggFine(gvFine);
+ auto glue = makeGlue(ggFine, ggCoarse);
+
+ // a quadratic function on the grid
+ auto f = [](const auto& pos){ return pos.two_norm2(); };
+
+ // make bases
+ P2 basisCoarse(gvCoarse), basisFine(gvFine);
SolutionVector solCoarse, solFine;
- interpolate(*basisCoarse, solCoarse, f);
- interpolate(*basisFine, solFine, f);
- auto gfCoarse = makeDiscreteGlobalBasisFunction<R>(*basisCoarse, solCoarse);
- auto gfFine = makeDiscreteGlobalBasisFunction<R>(*basisFine, solFine);
+ interpolate(basisCoarse, solCoarse, f);
+ interpolate(basisFine, solFine, f);
+ auto gfCoarse = makeDiscreteGlobalBasisFunction<R>(basisCoarse, solCoarse);
+ auto gfFine = makeDiscreteGlobalBasisFunction<R>(basisFine, solFine);
- // project fine discrete function onto coarse grid
- Dumux::FEGridGeometry<P2> ggCoarse(basisCoarse), ggFine(basisFine);
- const auto projector = Dumux::makeProjector(*basisFine, *basisCoarse, makeGlue(ggFine, ggCoarse));
+ // project fine discrete function onto coarse grid and convert back to grid function
+ const auto projector = Dumux::makeProjector(basisFine, basisCoarse, glue);
SolutionVector solCoarse2;
projector.project(solFine, solCoarse2);
- auto gfCoarse2 = makeDiscreteGlobalBasisFunction<R>(*basisCoarse, solCoarse);
+ auto gfCoarse2 = makeDiscreteGlobalBasisFunction<R>(basisCoarse, solCoarse);
// check that these functions are identical and that they exactly represent the analytical function
auto l2norm = Dumux::integrateL2Error(gvCoarse, gfCoarse, gfCoarse2, 3);
if (l2norm > 1e-14)
- DUNE_THROW(Dune::MathError, "L2 norm of projection and interpolation is non-zero!");
+ DUNE_THROW(Dune::MathError, "L2 norm of projection and interpolation is non-zero (" << l2norm << ")");
l2norm = Dumux::integrateL2Error(gvCoarse, makeAnalyticGridViewFunction(f, gvCoarse), gfCoarse2, 3);
if (l2norm > 1e-14)
- DUNE_THROW(Dune::MathError, "L2 norm of analytical solution and interpolation is non-zero!");
+ DUNE_THROW(Dune::MathError, "L2 norm of analytical solution and interpolation is non-zero (" << l2norm << ")");
+
+ // check the FV-style interface of integrateL2Error
+ Dumux::BoxFVGridGeometry<double, GridView> ggBoxCoarse(gvCoarse), ggBoxFine(gvFine);
+ auto p1BasisCoarse = Dumux::getFunctionSpaceBasis(ggBoxCoarse);
+ auto p1BasisFine = Dumux::getFunctionSpaceBasis(ggBoxFine);
+ auto f2 = [](const auto& pos){ return pos[0]; };
+ interpolate(p1BasisFine, solFine, f2);
+ interpolate(p1BasisCoarse, solCoarse, f2);
+
+ // project fine discrete function onto coarse grid
+ const auto projector2 = Dumux::makeProjector(p1BasisFine, p1BasisCoarse, glue);
+ projector2.project(solFine, solCoarse2);
+
+ l2norm = Dumux::integrateL2Error(ggBoxCoarse, solCoarse, solCoarse2, 3);
+ if (l2norm > 1e-14)
+ DUNE_THROW(Dune::MathError, "L2 norm of FV-style projected solutions is non-zero (" << l2norm << ")");
std::cout << "All tests passed!" << std::endl;
return 0;
--
GitLab