diff --git a/dune/gdt/local/numerical-fluxes/kinetic.hh b/dune/gdt/local/numerical-fluxes/kinetic.hh
index bbd0564c938fe127d38c40f153d1538acca81105..04a4daaa0d5ec1742a8685229413aa9972b6dfe7 100644
--- a/dune/gdt/local/numerical-fluxes/kinetic.hh
+++ b/dune/gdt/local/numerical-fluxes/kinetic.hh
@@ -13,8 +13,8 @@
#include <dune/xt/la/container.hh>
-#include <dune/gdt/momentmodels/basisfunctions.hh>
-#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
#include "interface.hh"
diff --git a/dune/gdt/operators/reconstruction/slopes.hh b/dune/gdt/operators/reconstruction/slopes.hh
index 3d884adfc3eb673ed9546dd74952d86816ba998c..15c461b94dbda56ee5a4c6e9cbf4ce710e726b4d 100644
--- a/dune/gdt/operators/reconstruction/slopes.hh
+++ b/dune/gdt/operators/reconstruction/slopes.hh
@@ -28,8 +28,8 @@
#include <dune/xt/common/fvector.hh>
#include <dune/xt/common/parallel/threadstorage.hh>
-#include <dune/gdt/momentmodels/basisfunctions/partial_moments.hh>
-#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
namespace Dune {
namespace GDT {
diff --git a/dune/gdt/test/entropic-coords-mn-discretization.hh b/dune/gdt/test/entropic-coords-mn-discretization.hh
index 753f612353129217c2f4d1f1f6b8cba2b727fa4f..1b8e4f3743f78a3314d09a5c59fb93b7d800b4ba 100644
--- a/dune/gdt/test/entropic-coords-mn-discretization.hh
+++ b/dune/gdt/test/entropic-coords-mn-discretization.hh
@@ -25,10 +25,10 @@
#include <dune/gdt/operators/generic.hh>
#include <dune/gdt/operators/reconstruction/linear_kinetic.hh>
#include <dune/gdt/interpolations/default.hh>
-#include <dune/gdt/momentmodels/hessianinverter.hh>
-#include <dune/gdt/momentmodels/entropyflux_kineticcoords.hh>
-#include <dune/gdt/momentmodels/entropyflux.hh>
-#include <dune/gdt/momentmodels/density_evaluations.hh>
+#include <dune/gdt/test/momentmodels/hessianinverter.hh>
+#include <dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/density_evaluations.hh>
#include <dune/gdt/local/numerical-fluxes/kinetic.hh>
#include <dune/gdt/local/operators/advection-fv.hh>
#include <dune/gdt/spaces/l2/finite-volume.hh>
diff --git a/dune/gdt/test/mn-discretization.hh b/dune/gdt/test/mn-discretization.hh
index 4ca09106572e66b906cec32bcd200666c5ab00d9..5df56da736506327e8149683f1585a959b9fb906 100644
--- a/dune/gdt/test/mn-discretization.hh
+++ b/dune/gdt/test/mn-discretization.hh
@@ -23,7 +23,7 @@
#include <dune/gdt/operators/advection-fv-entropybased.hh>
#include <dune/gdt/operators/advection-fv.hh>
#include <dune/gdt/interpolations/default.hh>
-#include <dune/gdt/momentmodels/entropysolver.hh>
+#include <dune/gdt/test/momentmodels/entropysolver.hh>
#include <dune/gdt/local/numerical-fluxes/kinetic.hh>
#include <dune/gdt/local/operators/advection-fv.hh>
#include <dune/gdt/spaces/l2/finite-volume.hh>
diff --git a/dune/gdt/momentmodels/basisfunctions.hh b/dune/gdt/test/momentmodels/basisfunctions.hh
similarity index 100%
rename from dune/gdt/momentmodels/basisfunctions.hh
rename to dune/gdt/test/momentmodels/basisfunctions.hh
diff --git a/dune/gdt/momentmodels/basisfunctions/hatfunctions.hh b/dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh
similarity index 99%
rename from dune/gdt/momentmodels/basisfunctions/hatfunctions.hh
rename to dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh
index b9aebac16b6922dda787d355301cf0b10d3baf57..c40825ef597b33ec4c63e2398755db27b4a8c0e1 100644
--- a/dune/gdt/momentmodels/basisfunctions/hatfunctions.hh
+++ b/dune/gdt/test/momentmodels/basisfunctions/hatfunctions.hh
@@ -16,7 +16,7 @@
#include <dune/xt/common/fmatrix.hh>
-#include <dune/gdt/momentmodels/triangulation.hh>
+#include <dune/gdt/test/momentmodels/triangulation.hh>
#include "interface.hh"
diff --git a/dune/gdt/momentmodels/basisfunctions/interface.hh b/dune/gdt/test/momentmodels/basisfunctions/interface.hh
similarity index 99%
rename from dune/gdt/momentmodels/basisfunctions/interface.hh
rename to dune/gdt/test/momentmodels/basisfunctions/interface.hh
index 6d959c4f7304a1bbbb9d29ba9a5f1d499158d56c..70564ac742203d7aa5eadefc63eb48fef9140f4c 100644
--- a/dune/gdt/momentmodels/basisfunctions/interface.hh
+++ b/dune/gdt/test/momentmodels/basisfunctions/interface.hh
@@ -22,7 +22,7 @@
#include <dune/xt/data/spherical_quadratures.hh>
#include <dune/gdt/discretefunction/default.hh>
-#include <dune/gdt/momentmodels/triangulation.hh>
+#include <dune/gdt/test/momentmodels/triangulation.hh>
namespace Dune {
namespace GDT {
diff --git a/dune/gdt/momentmodels/basisfunctions/legendre.hh b/dune/gdt/test/momentmodels/basisfunctions/legendre.hh
similarity index 100%
rename from dune/gdt/momentmodels/basisfunctions/legendre.hh
rename to dune/gdt/test/momentmodels/basisfunctions/legendre.hh
diff --git a/dune/gdt/momentmodels/basisfunctions/partial_moments.hh b/dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh
similarity index 100%
rename from dune/gdt/momentmodels/basisfunctions/partial_moments.hh
rename to dune/gdt/test/momentmodels/basisfunctions/partial_moments.hh
diff --git a/dune/gdt/momentmodels/basisfunctions/spherical_harmonics.hh b/dune/gdt/test/momentmodels/basisfunctions/spherical_harmonics.hh
similarity index 100%
rename from dune/gdt/momentmodels/basisfunctions/spherical_harmonics.hh
rename to dune/gdt/test/momentmodels/basisfunctions/spherical_harmonics.hh
diff --git a/dune/gdt/momentmodels/density_evaluations.hh b/dune/gdt/test/momentmodels/density_evaluations.hh
similarity index 99%
rename from dune/gdt/momentmodels/density_evaluations.hh
rename to dune/gdt/test/momentmodels/density_evaluations.hh
index 113589e67ad4687b819ff6c9772a8512758c2fe3..c3827d7a29d3a154cab65826f04b40cc6c7d84d1 100644
--- a/dune/gdt/momentmodels/density_evaluations.hh
+++ b/dune/gdt/test/momentmodels/density_evaluations.hh
@@ -17,7 +17,7 @@
#include <dune/xt/common/parameter.hh>
#include <dune/gdt/discretefunction/default.hh>
-#include <dune/gdt/momentmodels/entropyflux_kineticcoords.hh>
+#include <dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh>
#include <dune/gdt/operators/interfaces.hh>
#include <dune/gdt/type_traits.hh>
diff --git a/dune/gdt/momentmodels/entropyflux.hh b/dune/gdt/test/momentmodels/entropyflux.hh
similarity index 99%
rename from dune/gdt/momentmodels/entropyflux.hh
rename to dune/gdt/test/momentmodels/entropyflux.hh
index d81860778b65daf9d1f51bd3f9788b4450443542..f9658b412b6d5df45d71d87906c04773f58c3ca6 100644
--- a/dune/gdt/momentmodels/entropyflux.hh
+++ b/dune/gdt/test/momentmodels/entropyflux.hh
@@ -21,8 +21,8 @@
#include <dune/xt/functions/interfaces/flux-function.hh>
-#include <dune/gdt/momentmodels/basisfunctions.hh>
-#include <dune/gdt/momentmodels/entropyflux_implementations.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions.hh>
+#include <dune/gdt/test/momentmodels/entropyflux_implementations.hh>
namespace Dune {
namespace GDT {
diff --git a/dune/gdt/momentmodels/entropyflux_implementations.hh b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
similarity index 99%
rename from dune/gdt/momentmodels/entropyflux_implementations.hh
rename to dune/gdt/test/momentmodels/entropyflux_implementations.hh
index b964d10d57ba18ba6c652d429a8ef318421f55c1..974a036fb3df558080c3f9820035afe33e04a591 100644
--- a/dune/gdt/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
@@ -40,7 +40,7 @@
#include <dune/xt/functions/interfaces/function.hh>
-#include <dune/gdt/momentmodels/basisfunctions.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions.hh>
#include <dune/gdt/type_traits.hh>
#include "config.h"
diff --git a/dune/gdt/momentmodels/entropyflux_kineticcoords.hh b/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
similarity index 98%
rename from dune/gdt/momentmodels/entropyflux_kineticcoords.hh
rename to dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
index 6db77042bce69a9e9ffa4af477c261c4cfd0bbef..02852a3e628dfd1f12929ed2aa2b118dde741299 100644
--- a/dune/gdt/momentmodels/entropyflux_kineticcoords.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
@@ -18,9 +18,9 @@
#include <dune/xt/functions/interfaces/flux-function.hh>
-#include <dune/gdt/momentmodels/basisfunctions.hh>
-#include <dune/gdt/momentmodels/entropyflux_implementations.hh>
-#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions.hh>
+#include <dune/gdt/test/momentmodels/entropyflux_implementations.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
namespace Dune {
namespace GDT {
diff --git a/dune/gdt/momentmodels/entropysolver.hh b/dune/gdt/test/momentmodels/entropysolver.hh
similarity index 99%
rename from dune/gdt/momentmodels/entropysolver.hh
rename to dune/gdt/test/momentmodels/entropysolver.hh
index 1c16194b2a140466858bff422eb767314fe68287..c3b34d0e68242d9e719b04f9270f24624257bc81 100644
--- a/dune/gdt/momentmodels/entropysolver.hh
+++ b/dune/gdt/test/momentmodels/entropysolver.hh
@@ -16,7 +16,7 @@
#include <dune/xt/common/parameter.hh>
#include <dune/gdt/discretefunction/default.hh>
-#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
#include <dune/gdt/operators/interfaces.hh>
#include <dune/gdt/type_traits.hh>
diff --git a/dune/gdt/momentmodels/hessianinverter.hh b/dune/gdt/test/momentmodels/hessianinverter.hh
similarity index 99%
rename from dune/gdt/momentmodels/hessianinverter.hh
rename to dune/gdt/test/momentmodels/hessianinverter.hh
index c37eeb04122afbccb33b7473366c2265e0c669ca..b3b6fee38fe721a013954406e14eb1b9d1760290 100644
--- a/dune/gdt/momentmodels/hessianinverter.hh
+++ b/dune/gdt/test/momentmodels/hessianinverter.hh
@@ -16,7 +16,7 @@
#include <dune/xt/common/parameter.hh>
#include <dune/gdt/discretefunction/default.hh>
-#include <dune/gdt/momentmodels/entropyflux_kineticcoords.hh>
+#include <dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh>
#include <dune/gdt/operators/interfaces.hh>
#include <dune/gdt/type_traits.hh>
diff --git a/dune/gdt/test/momentmodels/kinetictransport/base.hh b/dune/gdt/test/momentmodels/kinetictransport/base.hh
index 9367a38d563c198105b1be324140392b3d991b1a..0a5239f6f759846944b410034fe089d61c0426a8 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/base.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/base.hh
@@ -17,8 +17,8 @@
#include <dune/xt/la/solver.hh>
-#include <dune/gdt/momentmodels/basisfunctions.hh>
-#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
#include <dune/gdt/test/momentmodels/kineticequation.hh>
#include "../kineticequation.hh"
diff --git a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
index 260a26d5a6dc3b2625154a6b0fb272edb7aabf53..f16229a02abd021ac0eda054f83b60b15e16b604 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
@@ -13,7 +13,7 @@
#include <dune/grid/yaspgrid.hh>
-#include <dune/gdt/momentmodels/basisfunctions.hh>
+#include <dune/gdt/test/momentmodels/basisfunctions.hh>
#include <dune/gdt/spaces/l2/finite-volume.hh>
#include <dune/gdt/spaces/l2/discontinuous-lagrange.hh>
#include <dune/gdt/tools/timestepper/interface.hh>
diff --git a/dune/gdt/test/momentmodels/moment-approximation.hh b/dune/gdt/test/momentmodels/moment-approximation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..967a1c7d9bc722348ce25df1f3141f7d4f2f0676
--- /dev/null
+++ b/dune/gdt/test/momentmodels/moment-approximation.hh
@@ -0,0 +1,427 @@
+// This file is part of the dune-gdt project:
+// https://github.com/dune-community/dune-gdt
+// Copyright 2010-2016 dune-gdt developers and contributors. All rights reserved.
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// Authors:
+// Tobias Leibner (2016)
+
+#ifndef DUNE_GDT_TEST_HYPERBOLIC_MOMENT_APPROXIMATION_HH
+#define DUNE_GDT_TEST_HYPERBOLIC_MOMENT_APPROXIMATION_HH
+
+#include <chrono>
+#include <fstream>
+
+#include <dune/xt/common/string.hh>
+
+#include <dune/xt/data/quadratures/fekete.hh>
+
+#include <dune/xt/grid/gridprovider.hh>
+
+#include <dune/xt/la/container.hh>
+
+#include <dune/gdt/type_traits.hh>
+#include <dune/gdt/test/momentmodels/entropyflux.hh>
+
+namespace Dune {
+namespace GDT {
+
+
+bool is_empty(const std::string filename)
+{
+ std::ifstream file(filename);
+ return !file.good() || file.peek() == std::ifstream::traits_type::eof();
+}
+
+template <class MomentBasisType, class DiscreteFunctionType>
+struct MomentApproximation
+{
+ using SpaceType = typename DiscreteFunctionType::SpaceType;
+ using GridViewType = typename SpaceType::GridViewType;
+ using GridType = typename GridViewType::Grid;
+ using DomainType = typename MomentBasisType::DomainType;
+ using RangeFieldType = typename DiscreteFunctionType::RangeFieldType;
+ using DynamicRangeType = typename MomentBasisType::DynamicRangeType;
+ using QuadraturesType = typename MomentBasisType::QuadraturesType;
+ static const size_t dimDomain = MomentBasisType::dimDomain;
+ static const size_t dimRange = MomentBasisType::dimRange;
+
+ // Chooses the quadrature that is used to calculate l1, l2, linf error and do the visualization
+ // In 1d, we use a Gauss-Lobatto-Quadrature of order 197 (100 quadrature points) on each interval. For the
+ // fullmoments, we fix the number of intervals to 50. Thus we have a different quadrature for each model, but the high
+ // order of the quadratures should ensure that the quadrature error is negligible compared to the model error.
+ // Choosing a fixed quadrature for all models would not allow to align the quadrature to the interval boundaries,
+ // which potentially causes much higher errors especially in linf norm. In 3d, all spherical triangles are refinements
+ // of the initial octants, so we can use the same quadrature (on the finest triangulation) for all models. Note that
+ // this quadrature is not used for the optimization problems in the entropy-based models. For the optimization
+ // problems, we use the following model-adapted quadratures: PMM_n and HFM_n (1d): A Gauss-Lobatto quadrature of order
+ // 15 on each interval. M_N (1d): A Gauss-Lobatto quadrature of order 197 on each of the intervals [-1,0] and [0,1].
+ // PMM_n and HFM_n (3d): A Fekete quadrature of order 9 on each spherical triangle (order 15 for PMM_32 and HFM_6).
+ // M_N (3d): A product quadrature of order 2N+8 on the octants of the sphere.
+ template <size_t domainDim = dimDomain, bool anything = true>
+ struct QuadratureHelper;
+
+ template <bool anything>
+ struct QuadratureHelper<3, anything>
+ {
+ static std::unique_ptr<SphericalTriangulation<RangeFieldType>> create_triangulation()
+ {
+ return std::make_unique<SphericalTriangulation<RangeFieldType>>(MomentBasisType::num_refinements);
+ }
+
+ static QuadraturesType get(const SphericalTriangulation<RangeFieldType>& triangulation,
+ const int quadrature_refinements,
+ const bool visualization,
+ const int /*additional_refs*/)
+ {
+ QuadratureRule<RangeFieldType, 2> reference_quadrature_rule =
+ XT::Data::FeketeQuadrature<RangeFieldType>::get(visualization ? 1 : 9);
+ const int additional_refs = quadrature_refinements - static_cast<int>(MomentBasisType::num_refinements);
+ if (additional_refs < 0)
+ DUNE_THROW(InvalidStateException,
+ "The number of refinements for the quadrature has to be greater or equal than the number of "
+ "refinements the basisfunctions use!");
+ return triangulation.quadrature_rules(static_cast<size_t>(additional_refs), reference_quadrature_rule);
+ }
+ };
+
+ template <bool anything>
+ struct QuadratureHelper<1, anything>
+ {
+ static std::unique_ptr<SphericalTriangulation<RangeFieldType>> create_triangulation()
+ {
+ return std::make_unique<SphericalTriangulation<RangeFieldType>>();
+ }
+
+ static QuadraturesType get(const SphericalTriangulation<RangeFieldType>& /*triangulation*/,
+ const int num_quadrature_intervals,
+ const bool /*visualization*/,
+ const int additional_refs)
+ {
+ return MomentBasisType::gauss_lobatto_quadratures(num_quadrature_intervals, 197, additional_refs);
+ }
+ };
+
+ // use block-wise solver for 3d partial moments, sparse solver for full matrix otherwise
+ template <bool partialmoments = is_3d_partial_moment_basis<MomentBasisType>::value, bool anything = true>
+ struct SolverHelper;
+
+ template <bool anything>
+ struct SolverHelper<true, anything>
+ {
+ static void solve(const MomentBasisType& basis_functions, const DynamicRangeType& u, DynamicRangeType& pn_coeffs)
+ {
+ const auto mass_matrix = basis_functions.block_mass_matrix();
+ FieldVector<double, 4> local_coeffs;
+ FieldVector<double, 4> local_u;
+ for (size_t jj = 0; jj < mass_matrix->num_rows / 4; ++jj) {
+ for (size_t kk = 0; kk < 4; ++kk)
+ local_u[kk] = u[4 * jj + kk];
+ mass_matrix->block(jj).solve(local_coeffs, local_u);
+ for (size_t kk = 0; kk < 4; ++kk)
+ pn_coeffs[4 * jj + kk] = local_coeffs[kk];
+ }
+ }
+ };
+
+ template <bool anything>
+ struct SolverHelper<false, anything>
+ {
+ static void solve(const MomentBasisType& basis_functions, const DynamicRangeType& u, DynamicRangeType& pn_coeffs)
+ {
+ const auto mass_matrix = basis_functions.mass_matrix();
+ XT::LA::IstlRowMajorSparseMatrix<double> sparse_mass_matrix(mass_matrix, true);
+ XT::LA::IstlDenseVector<double> pn_coeffs_istl(sparse_mass_matrix.rows());
+ XT::LA::IstlDenseVector<double> u_istl(pn_coeffs_istl);
+ for (size_t ii = 0; ii < u_istl.size(); ++ii)
+ u_istl.set_entry(ii, u[ii]);
+ XT::LA::solve(sparse_mass_matrix, u_istl, pn_coeffs_istl);
+ for (size_t ii = 0; ii < pn_coeffs_istl.size(); ++ii)
+ pn_coeffs[ii] = pn_coeffs_istl.get_entry(ii);
+ }
+ };
+
+ static void run(const int quad_refinements, std::string testcase, std::string filename = "")
+ {
+ //******************* create grid and basisfunctions ***************************************
+ auto grid_config = Dune::XT::Grid::cube_gridprovider_default_config();
+ const auto grid_provider =
+ Dune::XT::Grid::CubeGridProviderFactory<GridType>::create(grid_config, MPIHelper::getCommunicator());
+ const auto grid_view = grid_provider.leaf_view();
+
+ // ***************** choose testcase *********************
+ RangeFieldType tau = 1e-9;
+ RangeFieldType additional_quad_refinements =
+ (dimDomain == 1 && quad_refinements < 50) ? std::ceil(std::log2(50. / quad_refinements)) : 0.;
+ std::function<RangeFieldType(const DomainType&, const bool)> psi;
+ if (testcase == "Heaviside") {
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ psi = [](const DomainType& v, const bool left) {
+ return v[0] < 0 || (XT::Common::is_zero(v[0]) && left) ? RangeFieldType(5e-9) : RangeFieldType(1);
+ };
+ } else if (testcase == "HeavisideUnfitted") {
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ psi = [](const DomainType& v, const bool /*left*/) {
+ return v[0] < 1. / M_PI ? RangeFieldType(5e-9) : RangeFieldType(1);
+ };
+ } else if (testcase == "Gauss1d") {
+ tau = 1e-12;
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const RangeFieldType sigma = 0.5;
+ const RangeFieldType mu = 0.;
+ psi = [sigma, mu](const DomainType& v, const bool) {
+ return 1. / std::sqrt(2. * M_PI * std::pow(sigma, 2))
+ * std::exp(std::pow(v[0] - mu, 2) / (-2 * std::pow(sigma, 2)));
+ };
+ } else if (testcase == "CrossingBeams1dSmooth") {
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const RangeFieldType factor = 1e5;
+ const RangeFieldType norm = std::sqrt(M_PI / factor);
+ psi = [factor, norm](const DomainType& v, const bool) {
+ return std::max((std::exp(-factor * std::pow(v[0] - 1, 2)) + std::exp(-factor * std::pow(v[0] + 1, 2))) / norm,
+ 5e-9);
+ };
+ } else if (testcase == "CrossingBeams1dSmooth_2") {
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const RangeFieldType factor = 1e5;
+ const RangeFieldType norm = std::sqrt(M_PI / factor);
+ psi = [factor, norm](const DomainType& v, const bool) {
+ return std::max(
+ (std::exp(-factor * std::pow(v[0] - 0.5, 2)) + std::exp(-factor * std::pow(v[0] + 1, 2))) / norm, 5e-9);
+ };
+ } else if (testcase == "CrossingBeams1dDiscontinuous") {
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ psi = [](const DomainType& v, const bool) {
+ RangeFieldType height = 100;
+ if (std::abs(v[0] + 1) < 1. / (2 * height) || std::abs(v[0] - 1) < 1. / (2 * height))
+ return height;
+ else
+ return 5e-9;
+ };
+ } else if (testcase == "CrossingBeams1dDiscontinuous_2") {
+ if (dimDomain != 1)
+ DUNE_THROW(InvalidStateException,
+ "This is a 1-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ psi = [](const DomainType& v, const bool) {
+ RangeFieldType height = 100;
+ if (std::abs(v[0] + 1) < 1. / (2 * height) || std::abs(v[0] - 0.5) < 1. / (2 * height))
+ return height;
+ else
+ return 5e-9;
+ };
+ } else if (testcase == "GaussOnSphere") {
+ tau = 1e-12;
+ if (dimDomain != 3)
+ DUNE_THROW(InvalidStateException,
+ "This is a 3-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const RangeFieldType sigma = 0.5;
+ const RangeFieldType mu = 1.;
+ psi = [sigma, mu](const DomainType& v, const bool) {
+ return 1. / (2. * M_PI * std::pow(sigma, 2) * (1. - std::exp(-2 / std::pow(sigma, 2))))
+ * std::exp((std::pow(v[0] - mu, 2) + std::pow(v[1], 2) + std::pow(v[2], 2)) / (-2 * std::pow(sigma, 2)));
+ };
+ } else if (testcase == "SquareOnSphere") {
+ if (dimDomain != 3)
+ DUNE_THROW(InvalidStateException,
+ "This is a 3-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ psi = [](const DomainType& v, const bool) {
+ if (v[0] > 0 && std::abs(v[1]) < 0.5 && std::abs(v[2]) < 0.5)
+ return RangeFieldType(1);
+ else
+ return RangeFieldType(1e-8 / (4 * M_PI));
+ };
+ } else if (testcase == "CrossingBeams3d") {
+ if (dimDomain != 3)
+ DUNE_THROW(InvalidStateException,
+ "This is a 3-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const RangeFieldType factor = 1e2;
+ const RangeFieldType norm = M_PI / factor;
+ const FieldVector<RangeFieldType, 3> center1{{1., 0., 0.}};
+ const FieldVector<RangeFieldType, 3> center2{{0., 1., 0.}};
+ psi = [factor, norm, center1, center2](const DomainType& v, const bool) {
+ return std::max((std::exp(-factor * (v - center1).two_norm2()) + std::exp(-factor * (v - center2).two_norm2()))
+ / norm,
+ 1e-8 / (4 * M_PI));
+ };
+ } else if (testcase == "CrossingBeams3d_2") {
+ if (dimDomain != 3)
+ DUNE_THROW(InvalidStateException,
+ "This is a 3-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const RangeFieldType factor = 1e2;
+ const RangeFieldType norm = M_PI / factor;
+ const FieldVector<RangeFieldType, 3> center1{{std::sqrt(1. / 3.), std::sqrt(1. / 3.), std::sqrt(1. / 3.)}};
+ const FieldVector<RangeFieldType, 3> center2{
+ {std::sqrt(1. / (2 * M_PI)), -std::sqrt(0.5 - 1. / (4 * M_PI)), std::sqrt(0.5 - 1. / (4 * M_PI))}};
+ psi = [factor, norm, center1, center2](const DomainType& v, const bool) {
+ return std::max((std::exp(-factor * (v - center1).two_norm2()) + std::exp(-factor * (v - center2).two_norm2()))
+ / norm,
+ 1e-8 / (4 * M_PI));
+ };
+ } else if (testcase == "CrossingBeams3dDiscontinuous") {
+ if (dimDomain != 3)
+ DUNE_THROW(InvalidStateException,
+ "This is a 3-dimensional test, but the basis functions are " + XT::Common::to_string(dimDomain)
+ + "-dimensional!");
+ const FieldVector<RangeFieldType, 3> center1{{std::sqrt(1. / 3.), std::sqrt(1. / 3.), std::sqrt(1. / 3.)}};
+ const FieldVector<RangeFieldType, 3> center2{
+ {std::sqrt(1. / (2 * M_PI)), -std::sqrt(0.5 - 1. / (4 * M_PI)), std::sqrt(0.5 - 1. / (4 * M_PI))}};
+ const RangeFieldType r = 0.1;
+ const RangeFieldType r_squared = std::pow(r, 2);
+ psi = [r_squared, center1, center2](const DomainType& v, const bool) {
+ if ((v - center1).two_norm2() < r_squared || (v - center2).two_norm2() < r_squared)
+ return 1 / (M_PI * r_squared);
+ else
+ return 1e-8 / (4 * M_PI);
+ };
+ } else {
+ DUNE_THROW(NotImplemented, "Unknown testcase " + testcase + "!");
+ }
+
+ //******************* choose quadrature that is used for visualization and error calculation *******************
+ auto begin = std::chrono::steady_clock::now();
+ const auto triangulation = QuadratureHelper<>::create_triangulation();
+ std::chrono::duration<double> time = std::chrono::steady_clock::now() - begin;
+ if (dimDomain == 3)
+ std::cout << "Creating triangulation took: " << time.count() << " s!" << std::endl;
+ begin = std::chrono::steady_clock::now();
+ const auto visualization_quadratures =
+ QuadratureHelper<>::get(*triangulation, quad_refinements, true, additional_quad_refinements);
+ const auto basis_quadratures =
+ QuadratureHelper<>::get(*triangulation, quad_refinements, false, additional_quad_refinements);
+ time = std::chrono::steady_clock::now() - begin;
+ std::cout << "Creating quadratures took: " << time.count() << " s!" << std::endl;
+ begin = std::chrono::steady_clock::now();
+ std::shared_ptr<const MomentBasisType> basis_functions =
+ std::make_shared<const MomentBasisType>(*triangulation, basis_quadratures);
+ time = std::chrono::steady_clock::now() - begin;
+ std::cout << "Creating basis functions took: " << time.count() << " s!" << std::endl;
+
+ begin = std::chrono::steady_clock::now();
+ const auto u = basis_functions->get_moment_vector(psi);
+ time = std::chrono::steady_clock::now() - begin;
+ std::cout << "Calculating u took: " << time.count() << " s!" << std::endl;
+ using MnFluxType = EntropyBasedFluxFunction<GridViewType, MomentBasisType>;
+ begin = std::chrono::steady_clock::now();
+ MnFluxType mn_flux(grid_view,
+ *basis_functions,
+ tau,
+ true,
+ 0.01,
+ 0.5,
+ 1e-3,
+ {0, 1e-8, 1e-6, 1e-4, 1e-3, 1e-2, 5e-2, 0.1, 0.5, 1},
+ 50,
+ 100);
+ std::cout << "Creating entropy flux took: "
+ << std::chrono::duration_cast<std::chrono::seconds>(std::chrono::steady_clock::now() - begin).count()
+ << "s!" << std::endl;
+ const auto mn_begin = std::chrono::steady_clock::now();
+ const auto mn_ret = mn_flux.get_alpha(u, true);
+ const auto mn_end = std::chrono::steady_clock::now();
+ const std::chrono::duration<double> mn_time = mn_end - mn_begin;
+ std::cout << "Mn solve took: " << mn_time.count() << " s!" << std::endl;
+
+ const auto alpha = mn_ret->first;
+
+ const auto entropy_string =
+ "_"
+ + std::string(MomentBasisType::entropy == EntropyType::MaxwellBoltzmann ? "MaxwellBoltzmann" : "BoseEinstein")
+ + "_";
+ const auto filename_mn = filename + entropy_string + basis_functions->mn_name();
+ const auto filename_pn = filename + "_" + basis_functions->pn_name();
+ std::ofstream mn_file(filename_mn + ".txt");
+ std::ofstream pn_file(filename_pn + ".txt");
+ const std::string mn_errors_filename = filename + entropy_string + basis_functions->short_id() + "m_errors.txt";
+ const std::string pn_errors_filename = filename + "_" + basis_functions->short_id() + "p_errors.txt";
+ std::ofstream mn_errors_file(mn_errors_filename, std::ios_base::app);
+ std::ofstream pn_errors_file(pn_errors_filename, std::ios_base::app);
+ if (is_empty(mn_errors_filename))
+ mn_errors_file << "n l1error l2error linferror time(s) r" << std::endl;
+ if (is_empty(pn_errors_filename))
+ pn_errors_file << "n l1error l2error linferror time(s)" << std::endl;
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ mn_file << "v" << dd << " ";
+ pn_file << "v" << dd << " ";
+ }
+ mn_file << "psi" << std::endl;
+ pn_file << "psi" << std::endl;
+ RangeFieldType l1norm(0), l2norm(0), linfnorm(0);
+ RangeFieldType l1error_mn(0), l2error_mn(0), linferror_mn(0);
+ RangeFieldType l1error_pn(0), l2error_pn(0), linferror_pn(0);
+ DynamicRangeType pn_coeffs(dimRange);
+ const auto pn_begin = std::chrono::steady_clock::now();
+ SolverHelper<>::solve(*basis_functions, u, pn_coeffs);
+ const auto pn_end = std::chrono::steady_clock::now();
+ const std::chrono::duration<double> pn_time = pn_end - pn_begin;
+ std::cout << "Pn solve took: " << pn_time.count() << " s!" << std::endl;
+ const auto quadrature = XT::Data::merged_quadrature(visualization_quadratures);
+ begin = std::chrono::steady_clock::now();
+ for (auto it = quadrature.begin(); it != quadrature.end(); ++it) {
+ const auto& quad_point = *it;
+ const auto v = quad_point.position();
+ const auto weight = quad_point.weight();
+ const auto basis = basis_functions->evaluate(v, it.first_index());
+ const auto psi_mn_prod = std::inner_product(basis.begin(), basis.end(), alpha.begin(), 0.);
+ const auto psi_mn = (MomentBasisType::entropy == EntropyType::MaxwellBoltzmann)
+ ? std::exp(psi_mn_prod)
+ : std::exp(psi_mn_prod) / (1 - std::exp(psi_mn_prod));
+ const auto psi_pn = pn_coeffs * basis;
+ for (size_t dd = 0; dd < dimDomain; ++dd) {
+ mn_file << XT::Common::to_string(v[dd], 15) << " ";
+ pn_file << XT::Common::to_string(v[dd], 15) << " ";
+ }
+ mn_file << XT::Common::to_string(psi_mn, 15) << std::endl;
+ pn_file << XT::Common::to_string(psi_pn, 15) << std::endl;
+ l1norm += std::abs(psi(v, basis_functions->is_negative(it))) * weight;
+ l2norm += std::pow(psi(v, basis_functions->is_negative(it)), 2) * weight;
+ linfnorm = std::max(std::abs(psi(v, basis_functions->is_negative(it))), linfnorm);
+ l1error_mn += std::abs(psi_mn - psi(v, basis_functions->is_negative(it))) * weight;
+ l2error_mn += std::pow(psi_mn - psi(v, basis_functions->is_negative(it)), 2) * weight;
+ linferror_mn = std::max(std::abs(psi_mn - psi(v, basis_functions->is_negative(it))), linferror_mn);
+ l1error_pn += std::abs(psi_pn - psi(v, basis_functions->is_negative(it))) * weight;
+ l2error_pn += std::pow(psi_pn - psi(v, basis_functions->is_negative(it)), 2) * weight;
+ linferror_pn = std::max(std::abs(psi_pn - psi(v, basis_functions->is_negative(it))), linferror_pn);
+ }
+ time = std::chrono::steady_clock::now() - begin;
+ std::cout << "Calculating errors took: " << time.count() << " s!" << std::endl;
+ mn_errors_file << dimRange << " " << XT::Common::to_string(l1error_mn, 15) << " "
+ << XT::Common::to_string(l2error_mn, 15) << " " << XT::Common::to_string(linferror_mn, 15) << " "
+ << mn_time.count() << " " << mn_ret->second.second << std::endl;
+ pn_errors_file << dimRange << " " << XT::Common::to_string(l1error_pn, 15) << " "
+ << XT::Common::to_string(l2error_pn, 15) << " " << XT::Common::to_string(linferror_pn, 15) << " "
+ << pn_time.count() << std::endl;
+ std::cout << basis_functions->mn_name() << " done!" << std::endl;
+ }
+};
+
+template <class MomentBasisType, class DiscreteFunctionType>
+const size_t MomentApproximation<MomentBasisType, DiscreteFunctionType>::dimDomain;
+template <class MomentBasisType, class DiscreteFunctionType>
+const size_t MomentApproximation<MomentBasisType, DiscreteFunctionType>::dimRange;
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_TEST_HYPERBOLIC_MOMENT_APPROXIMATION_HH
diff --git a/dune/gdt/momentmodels/triangulation.hh b/dune/gdt/test/momentmodels/triangulation.hh
similarity index 100%
rename from dune/gdt/momentmodels/triangulation.hh
rename to dune/gdt/test/momentmodels/triangulation.hh
diff --git a/dune/gdt/type_traits.hh b/dune/gdt/type_traits.hh
index 220064eaf1100d5bd290ec9015e7f0ba8e03221d..cb8e7abf3d07528647c47230a2221323194fd9de 100644
--- a/dune/gdt/type_traits.hh
+++ b/dune/gdt/type_traits.hh
@@ -80,7 +80,7 @@ class LocalFiniteElementFamilyInterface;
template <class GV, size_t r, size_t rC, class R>
class SpaceInterface;
-// from #include <dune/gdt/momentmodels/basisfunctions.hh>
+// from #include <dune/gdt/test/momentmodels/basisfunctions.hh>
enum class EntropyType;
template <class DomainFieldType,
@@ -173,7 +173,7 @@ struct is_space<S, true> : public std::is_base_of<SpaceInterface<typename S::GV,
{};
-// from #include <dune/gdt/momentmodels/basisfunctions.hh>
+// from #include <dune/gdt/test/momentmodels/basisfunctions.hh>
template <class T>
struct is_hatfunction_basis : public std::false_type
{};