diff --git a/dune/gdt/test/momentmodels/density_evaluations.hh b/dune/gdt/test/momentmodels/density_evaluations.hh
index 6258f9996defd9ae047ceeb49ebb9dd728da0cec..fa8ed97f5b67cd2b93a2afd4f4b54fb3795d37b6 100644
--- a/dune/gdt/test/momentmodels/density_evaluations.hh
+++ b/dune/gdt/test/momentmodels/density_evaluations.hh
@@ -140,7 +140,9 @@ public:
     , space_(space)
     , boundary_distribution_(boundary_distribution)
     , min_acceptable_density_(min_acceptable_density)
-  {}
+  {
+    analytical_flux_.prepare_storage(space_.grid_view());
+  }
 
   bool linear() const override final
   {
@@ -159,13 +161,6 @@ public:
 
   void apply(const VectorType& alpha, VectorType& range, const XT::Common::Parameter& param = {}) const override final
   {
-    const auto num_entities = space_.grid_view().size(0);
-    analytical_flux_.exp_evaluations().resize(num_entities);
-    if (EntropyFluxType::entropy != EntropyType::MaxwellBoltzmann) {
-      analytical_flux_.eta_ast_prime_evaluations().resize(num_entities);
-      analytical_flux_.eta_ast_twoprime_evaluations().resize(num_entities);
-    }
-    analytical_flux_.boundary_distribution_evaluations().resize(num_entities);
     LocalDensityEvaluatorType local_density_evaluator(
         space_, alpha, range, analytical_flux_, boundary_distribution_, min_acceptable_density_, param);
     auto walker = XT::Grid::Walker<typename SpaceType::GridViewType>(space_.grid_view());
diff --git a/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh b/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
index f6b4529cb7562303a4044712684c4ba3040c8194..8710befd6e0c58e6a54b5a373aa432e91b1f018b 100644
--- a/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
+++ b/dune/gdt/test/momentmodels/entropic-coords-mn-discretization.hh
@@ -31,6 +31,7 @@
 #include <dune/gdt/spaces/l2/finite-volume.hh>
 #include <dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh>
 #include <dune/gdt/test/momentmodels/entropyflux.hh>
+#include <dune/gdt/test/momentmodels/entropysolver.hh>
 #include <dune/gdt/test/momentmodels/hessianinverter.hh>
 #include <dune/gdt/test/momentmodels/density_evaluations.hh>
 #include <dune/gdt/tools/timestepper/adaptive-rungekutta-kinetic.hh>
@@ -128,12 +129,19 @@ struct HyperbolicEntropicCoordsMnDiscretization
 
     // ***************** project initial values to discrete function *********************
     // create a discrete function for the solution
-    DiscreteFunctionType u(fv_space, "solution");
-    DiscreteFunctionType alpha(fv_space, "solution");
+    DiscreteFunctionType u(fv_space, "u_initial");
+    DiscreteFunctionType alpha(fv_space, "alpha_initial");
     // project initial values
     default_interpolation(*initial_values_u, u, grid_view);
 
     // convert initial values to alpha
+    // using EntropySolverType = EntropySolver<MomentBasis, SpaceType, MatrixType>;
+    // EntropySolverType entropy_solver(*entropy_flux,
+    //                                  fv_space,
+    //                                  problem.psi_vac() * basis_functions->unit_ball_volume() / 10,
+    //                                  filename);
+    // entropy_solver.apply(u.dofs().vector(), alpha.dofs().vector(), {});
+
     const auto u_local_func = u.local_discrete_function();
     const auto alpha_local_func = alpha.local_discrete_function();
     XT::Common::FieldVector<RangeFieldType, dimRange> u_local;
@@ -364,8 +372,8 @@ struct HyperbolicEntropicCoordsMnTest
   void run()
   {
     auto norms = HyperbolicEntropicCoordsMnDiscretization<TestCaseType>::run(
-                     100,
-                     0,
+                     10,
+                     -1,
                      TestCaseType::quad_order,
                      TestCaseType::quad_refinements,
                      "",
diff --git a/dune/gdt/test/momentmodels/entropyflux_implementations.hh b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
index f31792700fda6d2969f5447168dc036826b060d4..c4d15804242013d8644dec0b1a589a3d4fa56669 100644
--- a/dune/gdt/test/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_implementations.hh
@@ -475,7 +475,6 @@ public:
   // stores evaluations of exp(alpha^T b(v_i)) for all quadrature points v_i
   void store_exp_evaluations(QuadratureWeightsType& exp_evaluations, const DomainType& alpha) const
   {
-    exp_evaluations.resize(quad_points_.size());
     this->calculate_scalar_products(alpha, M_, exp_evaluations);
     this->apply_exponential(exp_evaluations);
   }
@@ -498,7 +497,6 @@ public:
       QuadratureWeightsType& boundary_distribution_evaluations,
       const std::function<RangeFieldType(const FluxDomainType&)>& boundary_distribution) const
   {
-    boundary_distribution_evaluations.resize(quad_points_.size());
     for (size_t ll = 0; ll < quad_points_.size(); ++ll)
       boundary_distribution_evaluations[ll] = boundary_distribution(quad_points_[ll]);
   }
@@ -1956,12 +1954,9 @@ public:
       QuadratureWeightsType& boundary_distribution_evaluations,
       const std::function<RangeFieldType(const FluxDomainType&)>& boundary_distribution) const
   {
-    boundary_distribution_evaluations.resize(num_blocks);
-    for (size_t jj = 0; jj < num_blocks; ++jj) {
-      boundary_distribution_evaluations[jj].resize(quad_points_[jj].size());
+    for (size_t jj = 0; jj < num_blocks; ++jj)
       for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
         boundary_distribution_evaluations[jj][ll] = boundary_distribution(quad_points_[jj][ll]);
-    }
   }
 
 
@@ -2130,11 +2125,8 @@ public:
                                  const BasisValuesMatrixType& M,
                                  QuadratureWeightsType& scalar_products) const
   {
-    scalar_products.resize(num_blocks);
-    for (size_t jj = 0; jj < num_blocks; ++jj) {
-      scalar_products[jj].resize(quad_weights_[jj].size());
+    for (size_t jj = 0; jj < num_blocks; ++jj)
       calculate_scalar_products_block(jj, beta_in.block(jj), M[jj], scalar_products[jj]);
-    }
   }
 
   void apply_exponential(BlockQuadratureWeightsType& values) const
@@ -2390,7 +2382,6 @@ public:
       for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
         M_[jj][ll] = basis_functions_.evaluate_on_face(quad_points_[jj][ll], jj);
     } // jj
-
   } // constructor
 
 
@@ -2918,12 +2909,9 @@ public:
       QuadratureWeightsType& boundary_distribution_evaluations,
       const std::function<RangeFieldType(const FluxDomainType&)>& boundary_distribution) const
   {
-    boundary_distribution_evaluations.resize(num_faces_);
-    for (size_t jj = 0; jj < num_faces_; ++jj) {
-      boundary_distribution_evaluations[jj].resize(quad_points_[jj].size());
+    for (size_t jj = 0; jj < num_faces_; ++jj)
       for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
         boundary_distribution_evaluations[jj][ll] = boundary_distribution(quad_points_[jj][ll]);
-    }
   }
 
 
@@ -2954,7 +2942,7 @@ public:
                                                const FluxDomainType& n_ij,
                                                const size_t dd) const
   {
-    thread_local FieldVector<QuadratureWeightsType, 2> eta_ast_prime_vals;
+    thread_local std::array<QuadratureWeightsType, 2> eta_ast_prime_vals;
     eta_ast_prime_vals[0].resize(num_faces_);
     eta_ast_prime_vals[1].resize(num_faces_);
     for (size_t jj = 0; jj < num_faces_; ++jj) {
@@ -3093,14 +3081,22 @@ public:
     return work_vecs;
   }
 
-  bool all_positive(const QuadratureWeightsType& vals) const
+  void resize_quad_weights_type(QuadratureWeightsType& weights) const
   {
+    weights.resize(num_faces_);
     for (size_t jj = 0; jj < num_faces_; ++jj)
+      weights[jj].resize(quad_weights_[jj].size());
+  }
+
+  bool all_positive(const QuadratureWeightsType& vals) const
+  {
+    for (size_t jj = 0; jj < num_faces_; ++jj) {
       for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
         const auto val = vals[jj][ll];
         if (val < 0. || std::isinf(val) || std::isnan(val))
           return false;
-      }
+      } // ll
+    } // jj
     return true;
   }
 
@@ -3108,9 +3104,7 @@ public:
   {
     LocalVectorType local_alpha;
     const auto& faces = basis_functions_.triangulation().faces();
-    scalar_products.resize(num_faces_);
     for (size_t jj = 0; jj < num_faces_; ++jj) {
-      scalar_products[jj].resize(quad_weights_[jj].size());
       const auto& vertices = faces[jj]->vertices();
       for (size_t ii = 0; ii < 3; ++ii)
         local_alpha[ii] = alpha.get_entry(vertices[ii]->index());
@@ -4429,7 +4423,6 @@ public:
       const std::function<RangeFieldType(const FluxDomainType&)>& boundary_distribution) const
   {
     for (size_t jj = 0; jj < num_intervals; ++jj) {
-      boundary_distribution_evaluations[jj].resize(quad_points_[jj].size());
       for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll)
         boundary_distribution_evaluations[jj][ll] = boundary_distribution(quad_points_[jj][ll]);
     }
diff --git a/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh b/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
index 1da29acc543df216472bfd5379a7a63661cae12a..a4ecb753c5d90e6ac0f3cd48f0ad86592e5e52d2 100644
--- a/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
+++ b/dune/gdt/test/momentmodels/entropyflux_kineticcoords.hh
@@ -241,6 +241,32 @@ public:
     return exp_evaluations_;
   }
 
+  void prepare_storage(const GridViewType& grid_view)
+  {
+    const auto num_entities = grid_view.size(0);
+    exp_evaluations_.resize(num_entities);
+    if (entropy != EntropyType::MaxwellBoltzmann) {
+      eta_ast_prime_storage_.resize(num_entities);
+      eta_ast_twoprime_storage_.resize(num_entities);
+    }
+    boundary_distribution_evaluations_.resize(num_entities);
+    for (auto&& entity : Dune::elements(grid_view)) {
+      const auto entity_index = grid_view.indexSet().index(entity);
+      implementation_->resize_quad_weights_type(exp_evaluations_[entity_index]);
+      if (entropy != EntropyType::MaxwellBoltzmann) {
+        implementation_->resize_quad_weights_type(eta_ast_prime_storage_[entity_index]);
+        implementation_->resize_quad_weights_type(eta_ast_twoprime_storage_[entity_index]);
+      }
+      for (auto&& intersection : Dune::intersections(grid_view, entity)) {
+        if (intersection.boundary()) {
+          const auto intersection_index = intersection.indexInInside();
+          implementation_->resize_quad_weights_type(
+              boundary_distribution_evaluations_[entity_index][intersection_index / 2][intersection_index % 2]);
+        }
+      } // intersections
+    } // entities
+  } // void prepare_storage(...)
+
   std::vector<QuadratureWeightsType>& eta_ast_prime_evaluations()
   {
     return *eta_ast_prime_evaluations_;