From 5f8d149d85829aa21ea2465699835baa5a34f3ca Mon Sep 17 00:00:00 2001
From: Tobias Leibner <tobias.leibner@googlemail.com>
Date: Thu, 4 Jul 2019 13:32:38 +0200
Subject: [PATCH] [momentmodels] finish implementing other entropies

---
 .ci/shared                                    |   2 +-
 .../entropyflux_implementations.hh            | 133 ++++++++++--------
 .../test/entropic-coords-mn-discretization.hh |   3 -
 ...bolic__momentmodels__entropic_coords_mn.cc |   2 -
 ...perbolic__momentmodels__mn_boseeinstein.cc |  17 ++-
 .../kinetictransport/testcases.hh             |  69 ++++++---
 dune/gdt/test/pn-discretization.hh            |  15 +-
 7 files changed, 158 insertions(+), 83 deletions(-)

diff --git a/.ci/shared b/.ci/shared
index 1752cb361..7816a58fe 160000
--- a/.ci/shared
+++ b/.ci/shared
@@ -1 +1 @@
-Subproject commit 1752cb36168d9f675b1fa21fa874c8fe2024c200
+Subproject commit 7816a58fe88344f9e2dea1c028e9dfbf0aef800e
diff --git a/dune/gdt/momentmodels/entropyflux_implementations.hh b/dune/gdt/momentmodels/entropyflux_implementations.hh
index 35829b4f4..5c6575c0f 100644
--- a/dune/gdt/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/momentmodels/entropyflux_implementations.hh
@@ -2608,7 +2608,7 @@ public:
     const auto& faces = basis_functions_.triangulation().faces();
     for (size_t jj = 0; jj < num_faces_; ++jj) {
       const auto& vertices = faces[jj]->vertices();
-      local_u *= 0.;
+      std::fill(local_u.begin(), local_u.end(), 0.);
       for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
         const auto factor = eta_ast_prime_vals[jj][ll] * quad_weights_[jj][ll];
         for (size_t ii = 0; ii < 3; ++ii)
@@ -3040,18 +3040,18 @@ public:
 #if ENTROPY_FLUX_USE_1D_HATFUNCTIONS_SPECIALIZATION
 #  if ENTROPY_FLUX_1D_HATFUNCTIONS_USE_ANALYTICAL_INTEGRALS
 /**
- * Specialization of EntropyBasedFluxImplementation for 1D Hatfunctions (no change of basis, analytic integrals +
- * Taylor)
+ * Specialization of EntropyBasedFluxImplementation for 1D Hatfunctions with MaxwellBoltzmann entropy
+ * (no change of basis, analytic integrals + Taylor)
  */
-template <class D, class R, size_t dimRange>
-class EntropyBasedFluxImplementation<HatFunctionMomentBasis<D, 1, R, dimRange, 1>>
+template <class D, class R, size_t dimRange, EntropyType entropy>
+class EntropyBasedFluxImplementation<HatFunctionMomentBasis<D, 1, R, dimRange, 1, 1, entropy>>
   : public XT::Functions::FunctionInterface<dimRange, 1, dimRange, R>
 {
   using BaseType = typename XT::Functions::FunctionInterface<dimRange, 1, dimRange, R>;
   using ThisType = EntropyBasedFluxImplementation;
 
 public:
-  using MomentBasis = HatFunctionMomentBasis<D, 1, R, dimRange, 1>;
+  using MomentBasis = HatFunctionMomentBasis<D, 1, R, dimRange, 1, 1, entropy>;
   static const size_t dimFlux = MomentBasis::dimFlux;
   static const size_t basis_dimRange = dimRange;
   using typename BaseType::DomainFieldType;
@@ -3064,6 +3064,7 @@ public:
   using FluxDomainType = FieldVector<DomainFieldType, dimFlux>;
   using VectorType = XT::Common::FieldVector<RangeFieldType, basis_dimRange>;
   using AlphaReturnType = std::pair<VectorType, std::pair<DomainType, RangeFieldType>>;
+  static_assert(entropy == EntropyType::MaxwellBoltzmann, "Not implemented for other entropies");
 
   explicit EntropyBasedFluxImplementation(const MomentBasis& basis_functions,
                                           const RangeFieldType tau,
@@ -3090,26 +3091,17 @@ public:
     , max_taylor_order_(max_taylor_order)
   {}
 
-  static bool is_realizable(const DomainType& u)
-  {
-    for (const auto& u_i : u)
-      if (!(u_i > 0.) || std::isinf(u_i))
-        return false;
-    return true;
-  }
+
+  // ============================================================================================
+  // ============================= FunctionInterface methods ====================================
+  // ============================================================================================
+
 
   virtual int order(const XT::Common::Parameter& /*param*/) const override
   {
     return 1;
   }
 
-  VectorType get_isotropic_alpha(const DomainType& u) const
-  {
-    static const auto alpha_iso = basis_functions_.alpha_iso();
-    static const auto alpha_one = basis_functions_.alpha_one();
-    return alpha_iso + alpha_one * std::log(basis_functions_.density(u));
-  }
-
   virtual RangeReturnType evaluate(const DomainType& u,
                                    const XT::Common::Parameter& /*param*/ = {}) const override final
   {
@@ -3193,6 +3185,12 @@ public:
     calculate_J_Hinv(result[0], J_diag, J_subdiag, H_diag, H_subdiag);
   }
 
+
+  // ============================================================================================
+  // =============================== Kinetic fluxes =============================================
+  // ============================================================================================
+
+
   // calculate \sum_{i=1}^d < v_i m \psi > n_i, where n is the unit outer normal,
   // m is the basis function vector, \psi is the ansatz corresponding to u
   // and x, v, t are the space, velocity and time variable, respectively
@@ -3369,16 +3367,17 @@ public:
     return ret;
   } // DomainType evaluate_kinetic_flux(...)
 
+
+  // ============================================================================================
+  // ============ Evaluations of ansatz distribution, moments, hessian etc. =====================
+  // ============================================================================================
+
+
   std::unique_ptr<AlphaReturnType> get_alpha(const DomainType& u) const
   {
     return get_alpha(u, get_isotropic_alpha(u), true);
   }
 
-  DomainType get_u(const DomainType& alpha) const
-  {
-    return calculate_u(alpha, M_, M_);
-  }
-
   // returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
   std::unique_ptr<AlphaReturnType>
   get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const
@@ -3483,39 +3482,11 @@ public:
     return ret;
   } // ... get_alpha(...)
 
-  const MomentBasis& basis_functions() const
+  DomainType get_u(const DomainType& alpha) const
   {
-    return basis_functions_;
+    return calculate_u(alpha);
   }
 
-private:
-  RangeFieldType calculate_f(const VectorType& alpha_k, const VectorType& v) const
-  {
-    RangeFieldType ret(0);
-    for (size_t ii = 0; ii < dimRange - 1; ++ii) {
-      if (std::abs(alpha_k[ii + 1] - alpha_k[ii]) > taylor_tol_) {
-        ret += (v_points_[ii + 1] - v_points_[ii]) / (alpha_k[ii + 1] - alpha_k[ii])
-               * (std::exp(alpha_k[ii + 1]) - std::exp(alpha_k[ii]));
-      } else {
-        RangeFieldType update = 1.;
-        RangeFieldType result = 0.;
-        size_t ll = 1;
-        RangeFieldType base = alpha_k[ii + 1] - alpha_k[ii];
-        auto pow_frac = 1.;
-        while (ll <= max_taylor_order_ && XT::Common::FloatCmp::ne(update, 0.)) {
-          update = pow_frac;
-          result += update;
-          ++ll;
-          pow_frac *= base / ll;
-        }
-        assert(!(std::isinf(pow_frac) || std::isnan(pow_frac)));
-        ret += result * (v_points_[ii + 1] - v_points_[ii]) * std::exp(alpha_k[ii]);
-      }
-    } // ii
-    ret -= alpha_k * v;
-    return ret;
-  } // .. calculate_f(...)
-
   VectorType calculate_u(const VectorType& alpha_k) const
   {
     VectorType u(0);
@@ -3566,6 +3537,34 @@ private:
     return u;
   } // VectorType calculate_u(...)
 
+
+  RangeFieldType calculate_f(const VectorType& alpha_k, const VectorType& v) const
+  {
+    RangeFieldType ret(0);
+    for (size_t ii = 0; ii < dimRange - 1; ++ii) {
+      if (std::abs(alpha_k[ii + 1] - alpha_k[ii]) > taylor_tol_) {
+        ret += (v_points_[ii + 1] - v_points_[ii]) / (alpha_k[ii + 1] - alpha_k[ii])
+               * (std::exp(alpha_k[ii + 1]) - std::exp(alpha_k[ii]));
+      } else {
+        RangeFieldType update = 1.;
+        RangeFieldType result = 0.;
+        size_t ll = 1;
+        RangeFieldType base = alpha_k[ii + 1] - alpha_k[ii];
+        auto pow_frac = 1.;
+        while (ll <= max_taylor_order_ && XT::Common::FloatCmp::ne(update, 0.)) {
+          update = pow_frac;
+          result += update;
+          ++ll;
+          pow_frac *= base / ll;
+        }
+        assert(!(std::isinf(pow_frac) || std::isnan(pow_frac)));
+        ret += result * (v_points_[ii + 1] - v_points_[ii]) * std::exp(alpha_k[ii]);
+      }
+    } // ii
+    ret -= alpha_k * v;
+    return ret;
+  } // .. calculate_f(...)
+
   VectorType calculate_gradient(const VectorType& alpha_k, const VectorType& v) const
   {
     return calculate_u(alpha_k) - v;
@@ -3770,6 +3769,30 @@ private:
         std::swap(ret_ptr[jj * dimRange + ii], ret_ptr[ii * dimRange + jj]);
   } // void calculate_J_Hinv(...)
 
+
+  // ============================================================================================
+  // ================================== Helper functions ========================================
+  // ============================================================================================
+
+
+  const MomentBasis& basis_functions() const
+  {
+    return basis_functions_;
+  }
+
+  static bool is_realizable(const DomainType& u)
+  {
+    for (const auto& u_i : u)
+      if (!(u_i > 0.) || std::isinf(u_i))
+        return false;
+    return true;
+  }
+
+  VectorType get_isotropic_alpha(const DomainType& u) const
+  {
+    return basis_functions_.alpha_iso(basis_functions_.density(u));
+  }
+
   const MomentBasis& basis_functions_;
   const std::vector<RangeFieldType>& v_points_;
   const RangeFieldType tau_;
diff --git a/dune/gdt/test/entropic-coords-mn-discretization.hh b/dune/gdt/test/entropic-coords-mn-discretization.hh
index 71e162969..3870ca7a0 100644
--- a/dune/gdt/test/entropic-coords-mn-discretization.hh
+++ b/dune/gdt/test/entropic-coords-mn-discretization.hh
@@ -109,9 +109,6 @@ struct HyperbolicEntropicCoordsMnDiscretization
     auto flux = problem.flux();
     auto* entropy_flux = dynamic_cast<OldEntropyFluxType*>(flux.get());
     auto analytical_flux = std::make_unique<EntropyFluxType>(*entropy_flux);
-    // for Legendre polynomials and real spherical harmonics, the results are sensitive to the initial guess in the
-    // Newton algorithm. If the thread cache is enabled, the guess is different dependent on how many threads we are
-    // using, so for the tests we disable this cache to get reproducible results.
     const RangeFieldType CFL = problem.CFL();
 
     // calculate boundary values for alpha
diff --git a/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc b/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
index 5cc4bd054..217629bbe 100644
--- a/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
+++ b/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
@@ -23,12 +23,10 @@ using YaspGridTestCasesAll = testing::Types<
     // Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 20>, true, true>
     Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 7>, false, true>,
     Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 7>, true, true>,
-#if 0
     Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
     Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, true, true>,
     Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
     Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, true, true>,
-#endif
     Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
     Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, true, true>,
     Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, false, true>,
diff --git a/dune/gdt/test/hyperbolic__momentmodels__mn_boseeinstein.cc b/dune/gdt/test/hyperbolic__momentmodels__mn_boseeinstein.cc
index ff4e3f025..0823feb15 100644
--- a/dune/gdt/test/hyperbolic__momentmodels__mn_boseeinstein.cc
+++ b/dune/gdt/test/hyperbolic__momentmodels__mn_boseeinstein.cc
@@ -27,10 +27,25 @@ using YaspGridTestCasesAll = testing::Types<
     Dune::GDT::
         SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1, 1, entropy>, false>,
     Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1, 1, entropy>, true>,
+    Dune::GDT::
+        SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1, entropy>, false>,
+    Dune::GDT::
+        SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1, entropy>, true>,
     Dune::GDT::
         PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3, entropy>, false>,
     Dune::GDT::
-        PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3, entropy>, true>>;
+        PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3, entropy>, true>
+#if HAVE_CLP
+    ,
+    Dune::GDT::PointSourceMnTestCase<Yasp3,
+                                     Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2, 3, false, entropy>,
+                                     true>
+#endif
+#if HAVE_QHULL
+    ,
+    Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::PartialMomentBasis<double, 3, double, 0, 1, 3, 1, entropy>, true>
+#endif
+    >;
 
 TYPED_TEST_CASE(HyperbolicMnTest, YaspGridTestCasesAll);
 TYPED_TEST(HyperbolicMnTest, check)
diff --git a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
index bfadf3dc8..5fac9aca6 100644
--- a/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
+++ b/dune/gdt/test/momentmodels/kinetictransport/testcases.hh
@@ -276,13 +276,25 @@ struct SourceBeamMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8,
   static constexpr double tol = 1e-9;
 };
 
+template <bool reconstruct>
+struct SourceBeamMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1, EntropyType::BoseEinstein>,
+                                   reconstruct,
+                                   false>
+{
+  static constexpr double l1norm = reconstruct ? 0.33140398337940113 : 0.33140398338096477;
+  static constexpr double l2norm = reconstruct ? 0.45580284843165519 : 0.44483205570831974;
+  static constexpr double linfnorm = reconstruct ? 0.99172119511603896 : 0.98930804287194951;
+  static constexpr double tol = 1e-9;
+};
+
+
 template <bool reconstruct>
 struct SourceBeamMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, true>
 {
   static constexpr double l1norm = reconstruct ? 371.54588397717055 : 367.97988291905477;
   static constexpr double l2norm = reconstruct ? 236.4476851910448 : 235.54814675091959;
   static constexpr double linfnorm = reconstruct ? 210.63369526083264 : 208.81107020771216;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <bool reconstruct>
@@ -311,7 +323,7 @@ struct SourceBeamMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1, 1
   static constexpr double l1norm = reconstruct ? 254.20216502516391 : 270.74268779687191;
   static constexpr double l2norm = reconstruct ? 187.86036790841933 : 202.76054800096165;
   static constexpr double linfnorm = reconstruct ? 265.10790627160509 : 260.82089045524185;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
@@ -406,7 +418,7 @@ struct PlaneSourceMnExpectedResults<LegendreMomentBasis<double, double, 7>, reco
   static constexpr double l1norm = reconstruct ? 33.830651291425575 : 31.119878976551046;
   static constexpr double l2norm = reconstruct ? 24.726893737746675 : 23.385570207485049;
   static constexpr double linfnorm = reconstruct ? 19.113827924512311 : 19.113827924512311;
-  static constexpr double tol = 1e-7;
+  static constexpr double tol = 1e-5;
 };
 
 template <bool reconstruct>
@@ -421,10 +433,10 @@ struct PlaneSourceMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8,
 template <bool reconstruct>
 struct PlaneSourceMnExpectedResults<HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, reconstruct, true>
 {
-  static constexpr double l1norm = reconstruct ? 268.42786311776274 : 246.7429359648828;
-  static constexpr double l2norm = reconstruct ? 197.1506642519208 : 186.09403264481648;
+  static constexpr double l1norm = reconstruct ? 268.42768559247429 : 246.7429359648828;
+  static constexpr double l2norm = reconstruct ? 197.1506094198385 : 186.09403264481648;
   static constexpr double linfnorm = reconstruct ? 152.91062339609854 : 152.91062339609854;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <bool reconstruct>
@@ -442,7 +454,7 @@ struct PlaneSourceMnExpectedResults<PartialMomentBasis<double, 1, double, 8, 1,
   static constexpr double l1norm = reconstruct ? 144.19157186249112 : 135.86834797834712;
   static constexpr double l2norm = reconstruct ? 104.28938402311856 : 100.2359224660796;
   static constexpr double linfnorm = reconstruct ? 100.43185554232102 : 97.933985765677008;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 
@@ -622,13 +634,25 @@ struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, do
   static constexpr double tol = 1e-9;
 };
 
+template <bool reconstruct>
+struct PointSourceMnExpectedResults<
+    RealSphericalHarmonicsMomentBasis<double, double, 2, 3, false, EntropyType::BoseEinstein>,
+    reconstruct,
+    false>
+{
+  static constexpr double l1norm = reconstruct ? 1.0000013830443903 : 0.;
+  static constexpr double l2norm = reconstruct ? 2.6909504479323516 : 0.;
+  static constexpr double linfnorm = reconstruct ? 10.375951173911345 : 0.;
+  static constexpr double tol = 1e-9;
+};
+
 template <bool reconstruct>
 struct PointSourceMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, double, 2, 3>, reconstruct, true>
 {
   static constexpr double l1norm = reconstruct ? 1674.9008041695579 : 1585.7044225325101;
   static constexpr double l2norm = reconstruct ? 619.41343145125302 : 589.93299566257235;
   static constexpr double linfnorm = reconstruct ? 264.16080528868997 : 266.5453598444696;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <bool reconstruct>
@@ -661,10 +685,10 @@ struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0,
 template <bool reconstruct>
 struct PointSourceMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct, true>
 {
-  static constexpr double l1norm = reconstruct ? 743.43592672927934 : 683.47651349520368;
-  static constexpr double l2norm = reconstruct ? 279.11700895978396 : 258.45009663177717;
+  static constexpr double l1norm = reconstruct ? 818.73622981959204 : 781.1965079003387;
+  static constexpr double l2norm = reconstruct ? 301.48148670872598 : 289.21738528985526;
   static constexpr double linfnorm = reconstruct ? 125.7102296804655 : 125.71014355518233;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <bool reconstruct>
@@ -676,13 +700,24 @@ struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1,
   static constexpr double tol = 1e-9;
 };
 
+template <bool reconstruct>
+struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1, 3, 1, EntropyType::BoseEinstein>,
+                                    reconstruct,
+                                    false>
+{
+  static constexpr double l1norm = reconstruct ? 1.0000000829624796 : 0.;
+  static constexpr double l2norm = reconstruct ? 2.6983603000374528 : 0.;
+  static constexpr double linfnorm = reconstruct ? 10.391283146511036 : 0.;
+  static constexpr double tol = 1e-9;
+};
+
 template <bool reconstruct>
 struct PointSourceMnExpectedResults<PartialMomentBasis<double, 3, double, 0, 1, 3, 1>, reconstruct, true>
 {
   static constexpr double l1norm = reconstruct ? 1167.5985275432627 : 1126.5174600848904;
   static constexpr double l2norm = reconstruct ? 428.15220455351266 : 415.19451310103005;
   static constexpr double linfnorm = reconstruct ? 188.26590907577059 : 191.48910050886076;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
@@ -721,7 +756,7 @@ struct CheckerboardMnExpectedResults<RealSphericalHarmonicsMomentBasis<double, d
   static constexpr double l1norm = reconstruct ? 0. : 0.;
   static constexpr double l2norm = reconstruct ? 0. : 0.;
   static constexpr double linfnorm = reconstruct ? 0. : 0.;
-  static constexpr double tol = 1e-9;
+  static constexpr double tol = 1e-5;
 };
 
 template <bool reconstruct>
@@ -736,10 +771,10 @@ struct CheckerboardMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0
 template <bool reconstruct>
 struct CheckerboardMnExpectedResults<HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, reconstruct, true>
 {
-  static constexpr double l1norm = reconstruct ? 44760.517221844952 : 0.;
-  static constexpr double l2norm = reconstruct ? 2491.7035345029099 : 0.;
-  static constexpr double linfnorm = reconstruct ? 542.02861997383934 : 0.;
-  static constexpr double tol = 1e-9;
+  static constexpr double l1norm = reconstruct ? 42799.981949017187 : 0.;
+  static constexpr double l2norm = reconstruct ? 2318.887531040597 : 0.;
+  static constexpr double linfnorm = reconstruct ? 131.24293776745421 : 0.;
+  static constexpr double tol = 1e-5;
 };
 
 template <class GridImp, class MomentBasisImp, bool reconstruct, bool kinetic_scheme = false>
diff --git a/dune/gdt/test/pn-discretization.hh b/dune/gdt/test/pn-discretization.hh
index 3cd158c0f..806e08d57 100644
--- a/dune/gdt/test/pn-discretization.hh
+++ b/dune/gdt/test/pn-discretization.hh
@@ -135,10 +135,17 @@ template <class AnalyticalFluxType,
           class DomainFieldType,
           size_t dimDomain,
           class RangeFieldType,
-          size_t dimRange_or_refinements>
-struct EigenvectorWrapperChooser<
-    Dune::GDT::PartialMomentBasis<DomainFieldType, dimDomain, RangeFieldType, dimRange_or_refinements, 1, dimDomain>,
-    AnalyticalFluxType>
+          size_t dimRange_or_refinements,
+          Dune::GDT::EntropyType entropy>
+struct EigenvectorWrapperChooser<Dune::GDT::PartialMomentBasis<DomainFieldType,
+                                                               dimDomain,
+                                                               RangeFieldType,
+                                                               dimRange_or_refinements,
+                                                               1,
+                                                               dimDomain,
+                                                               1,
+                                                               entropy>,
+                                 AnalyticalFluxType>
 {
   using type = Dune::GDT::internal::BlockedEigenvectorWrapper<AnalyticalFluxType>;
 };
-- 
GitLab