[test] complete rewrite and modernize of InstationaryEocStudy

dune/gdt/test/instationary-eocstudies/base.hh

+// This file is part of the dune-gdt project:
+//   https://github.com/dune-community/dune-gdt
+// Copyright 2010-2018 dune-gdt developers and contributors. All rights reserved.
+// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
+//          with "runtime exception" (http://www.dune-project.org/license.html)
+// Authors:
+//   Felix Schindler (2016 - 2017)
+//   Rene Milk       (2016 - 2018)
+//   Tobias Leibner  (2016 - 2017)
+
+#ifndef DUNE_GDT_TEST_INSTATIONARY_EOCSTUDIES_BASE_HH
+#define DUNE_GDT_TEST_INSTATIONARY_EOCSTUDIES_BASE_HH
+
+#include <cmath>
+#include <functional>
+#include <memory>
+
+#include <dune/grid/io/file/dgfparser/dgfparser.hh>
+
+#include <dune/xt/common/convergence-study.hh>
+#include <dune/xt/common/fvector.hh>
+#include <dune/xt/common/test/gtest/gtest.h>
+#include <dune/xt/la/container.hh>
+#include <dune/xt/la/container/vector-array/list.hh>
+#include <dune/xt/la/type_traits.hh>
+#include <dune/xt/grid/gridprovider/provider.hh>
+#include <dune/xt/grid/type_traits.hh>
+#include <dune/xt/functions/lambda/function.hh>
+
+#include <dune/gdt/exceptions.hh>
+#include <dune/gdt/functionals/localizable-functional.hh>
+#include <dune/gdt/local/bilinear-forms/integrals.hh>
+#include <dune/gdt/local/functionals/integrals.hh>
+#include <dune/gdt/local/integrands/abs.hh>
+#include <dune/gdt/local/integrands/product.hh>
+#include <dune/gdt/operators/interfaces.hh>
+#include <dune/gdt/operators/localizable-bilinear-form.hh>
+#include <dune/gdt/prolongations.hh>
+#include <dune/gdt/spaces/bochner.hh>
+#include <dune/gdt/spaces/interface.hh>
+#include <dune/gdt/type_traits.hh>
+
+namespace Dune {
+namespace GDT {
+namespace Test {
+
+
+template <class GridView, size_t m_ = 1, XT::LA::Backends la = XT::LA::Backends::istl_sparse>
+class InstationaryEocStudy : public XT::Common::ConvergenceStudy, public ::testing::Test
+{
+  static_assert(XT::Grid::is_view<GridView>::value, "");
+
+protected:
+  using GV = GridView;
+  using G = typename GridView::Grid;
+  static const constexpr size_t m = m_;
+  using D = double;
+  static const constexpr size_t d = G::dimension;
+  using R = double;
+  using E = XT::Grid::extract_entity_t<GV>;
+  using DomainType = XT::Common::FieldVector<D, d>;
+  using RangeType = XT::Common::FieldVector<D, m>;
+  using GP = XT::Grid::GridProvider<G>;
+  using S = SpaceInterface<GV, m>;
+  using V = typename XT::LA::Container<R, la>::VectorType;
+  using DF = DiscreteFunction<V, GV, m>;
+  using BS = BochnerSpace<GV, m>;
+  using O = OperatorInterface<V, GV, m>;
+
+public:
+  InstationaryEocStudy(const double T_end,
+                       const size_t num_refinements = 3,
+                       const size_t num_additional_refinements_for_reference = 2,
+                       std::function<void(const DiscreteBochnerFunction<V, GV, m>&, const std::string&)> visualizer =
+                           [](const auto& /*solution*/, const auto& /*prefix*/) { /*no visualization by default*/ })
+    : T_end_(T_end)
+    , num_refinements_(num_refinements)
+    , num_additional_refinements_for_reference_(num_additional_refinements_for_reference)
+    , visualize_(visualizer)
+    , current_refinement_(std::numeric_limits<size_t>::max())
+    , current_data_()
+    , current_grid_(nullptr)
+    , current_space_(nullptr)
+    , reference_grid_(nullptr)
+    , reference_space_(nullptr)
+    , current_solution_on_reference_grid_(nullptr)
+    , reference_solution_on_reference_grid_(nullptr)
+  {
+  }
+
+  virtual size_t num_refinements() const override
+  {
+    return num_refinements_;
+  }
+
+  virtual std::vector<std::string> targets() const override
+  {
+    return {"h", "dt"};
+  }
+
+  virtual std::vector<std::string> norms() const override
+  {
+    return {"L_infty/L_1", "L_infty/L_2", "L_2/L_2"};
+  }
+
+  virtual std::vector<std::string> estimates() const override
+  {
+    return {};
+  }
+
+  virtual std::string discretization_info_title() const override
+  {
+    return " |grid| |   #DoFs";
+  }
+
+protected:
+  virtual GP make_initial_grid() = 0;
+
+  virtual std::unique_ptr<S> make_space(const GP& current_grid) = 0;
+
+  virtual double estimate_dt(const S& space) = 0;
+
+public:
+  virtual std::string discretization_info(const size_t refinement_level) override
+  {
+    auto& self = *this;
+    if (current_refinement_ != refinement_level) {
+      // clear the current state
+      current_grid_.reset();
+      current_space_.reset();
+      current_solution_on_reference_grid_.reset();
+      // compute on this refinement
+      current_grid_ = std::make_unique<GP>(make_initial_grid());
+      for (size_t ref = 0; ref < refinement_level; ++ref)
+        current_grid_->global_refine(DGFGridInfo<G>::refineStepsForHalf());
+      current_space_ = make_space(*current_grid_);
+      double grid_size = 0;
+      double grid_width = 0.;
+      for (auto&& grid_element : elements(current_space_->grid_view())) {
+        grid_size += 1;
+        grid_width = std::max(grid_width, XT::Grid::entity_diameter(grid_element));
+      }
+      current_refinement_ = refinement_level;
+      current_data_.clear();
+      current_data_["target"]["h"] = grid_width;
+      current_data_["target"]["dt"] = self.estimate_dt(*current_space_);
+      current_data_["info"]["num_grid_elements"] = grid_size;
+      current_data_["info"]["num_dofs"] = current_space_->mapper().size();
+    }
+    const auto lfill_nicely = [&](const auto& number, const auto& len) {
+      std::string ret;
+      using XT::Common::to_string;
+      if (std::log10(number) < len)
+        ret = this->lfill(to_string(number), len);
+      else {
+        std::stringstream ss;
+        ss << std::setprecision(2) << std::scientific << number;
+        ret = this->lfill(ss.str(), len);
+      }
+      return ret;
+    };
+    return lfill_nicely(current_data_["info"]["num_grid_elements"], 7) + " | "
+           + lfill_nicely(current_data_["info"]["num_dofs"], 7);
+  } // ... discretization_info(...)
+
+protected:
+  virtual std::unique_ptr<O> make_lhs_operator(const S& space) = 0;
+
+  virtual XT::LA::ListVectorArray<V> solve(const S& space, const double T_end, const double dt) = 0;
+
+  virtual void compute_reference_solution()
+  {
+    if (reference_solution_on_reference_grid_)
+      return;
+    reference_grid_ = std::make_unique<GP>(make_initial_grid());
+    for (size_t ref = 0; ref < num_refinements_ + num_additional_refinements_for_reference_; ++ref)
+      reference_grid_->global_refine(DGFGridInfo<G>::refineStepsForHalf());
+    reference_space_ = make_space(*reference_grid_);
+    reference_solution_on_reference_grid_ =
+        std::make_unique<XT::LA::ListVectorArray<V>>(solve(*reference_space_, T_end_, estimate_dt(*reference_space_)));
+    // visualize
+    const BS reference_bochner_space(*reference_space_,
+                                     time_points_from_vector_array(*reference_solution_on_reference_grid_));
+    visualize_(make_discrete_bochner_function(reference_bochner_space, *reference_solution_on_reference_grid_),
+               "reference_solution_on_refinement_"
+                   + XT::Common::to_string(num_refinements_ + num_additional_refinements_for_reference_));
+  } // ... compute_reference_solution(...)
+
+public:
+  virtual std::map<std::string, std::map<std::string, double>>
+  compute(const size_t refinement_level, const std::vector<std::string>& only_these) override
+  {
+    auto& self = *this;
+    if (current_refinement_ != refinement_level)
+      self.discretization_info(refinement_level);
+    DUNE_THROW_IF(!current_space_, InvalidStateException, "");
+    // compute current solution
+    const auto& current_space = *current_space_;
+    Timer timer;
+    auto solution_on_current_grid = solve(current_space, T_end_, current_data_["target"]["dt"]);
+    current_data_["info"]["time_to_solution"] = timer.elapsed();
+    // visualize
+    const BS current_bochner_space(current_space, time_points_from_vector_array(solution_on_current_grid));
+    visualize_(make_discrete_bochner_function(current_bochner_space, solution_on_current_grid),
+               "solution_on_refinement_" + XT::Common::to_string(refinement_level));
+    const auto coarse_solution = make_discrete_bochner_function(current_bochner_space, solution_on_current_grid);
+    // determine wether we need a reference solution
+    // compute statistics
+    std::vector<std::string> actual_norms = self.filter(self.norms(), only_these);
+    // - norms
+    if (!actual_norms.empty()) {
+      self.compute_reference_solution();
+      DUNE_THROW_IF(!reference_space_, InvalidStateException, "");
+      const auto& reference_space = *reference_space_;
+      DUNE_THROW_IF(!reference_solution_on_reference_grid_, InvalidStateException, "");
+      auto& u = *reference_solution_on_reference_grid_;
+      // prolong
+      const BS reference_bochner_space(reference_space,
+                                       time_points_from_vector_array(*reference_solution_on_reference_grid_));
+      current_solution_on_reference_grid_ = std::make_unique<XT::LA::ListVectorArray<V>>(
+          prolong<V>(coarse_solution, reference_bochner_space).dof_vectors());
+      auto& u_h = *current_solution_on_reference_grid_;
+      while (!actual_norms.empty()) {
+        const auto norm_id = actual_norms.back();
+        const auto components = XT::Common::tokenize(norm_id, "/");
+        actual_norms.pop_back();
+        // compute Bochner norm
+        DUNE_THROW_IF(components.size() != 2,
+                      XT::Common::Exceptions::wrong_input_given,
+                      "I do not know how to compute the following norm: " << norm_id);
+        const auto temporal_norm_id = components[0];
+        const auto spatial_norm_id = components[1];
+        // - spatial component
+        std::function<double(const DF&)> spatial_norm;
+        if (spatial_norm_id == "L_1") {
+          spatial_norm = [&](const DF& func) {
+            auto localizable_functional = make_localizable_functional(reference_space.grid_view(), func);
+            localizable_functional.append(LocalElementIntegralFunctional<E>(LocalElementAbsIntegrand<E>()));
+            localizable_functional.assemble();
+            return localizable_functional.result();
+          };
+        } else if (spatial_norm_id == "L_2") {
+          spatial_norm = [&](const DF& func) {
+            auto localizable_product = make_localizable_bilinear_form(reference_space.grid_view(), func, func);
+            localizable_product.append(LocalElementIntegralBilinearForm<E>(LocalElementProductIntegrand<E>()));
+            localizable_product.assemble();
+            return std::sqrt(localizable_product.result());
+          };
+        } else
+          DUNE_THROW(XT::Common::Exceptions::wrong_input_given,
+                     "I do not know how to compute the spatial norm '" << spatial_norm_id << "'!");
+        // - temporal component
+        if (temporal_norm_id == "L_infty") {
+          double result = 0;
+          for (size_t ii = 0; ii < u.length(); ++ii)
+            result = std::max(result,
+                              spatial_norm(make_discrete_function(reference_space, u[ii].vector() - u_h[ii].vector())));
+          current_data_["norm"][norm_id] = result;
+        } else if (temporal_norm_id == "L_2") {
+          const XT::Functions::LambdaFunction<1> spatial_norm_function(1, [&](const auto& time, const auto& /*param*/) {
+            const auto u_t = make_discrete_bochner_function(reference_bochner_space, u).evaluate(time);
+            const auto u_h_t = make_discrete_bochner_function(reference_bochner_space, u_h).evaluate(time);
+            return spatial_norm(make_discrete_function(reference_space, u_t - u_h_t));
+          });
+          auto temporal_grid_view = reference_bochner_space.temporal_space().grid_view();
+          using TE = XT::Grid::extract_entity_t<decltype(temporal_grid_view)>;
+          const auto& spatial_norm_as_temporal_grid_function = spatial_norm_function.template as_grid_function<TE>();
+          auto localizable_temporal_product = make_localizable_bilinear_form(
+              temporal_grid_view, spatial_norm_as_temporal_grid_function, spatial_norm_as_temporal_grid_function);
+          localizable_temporal_product.append(LocalElementIntegralBilinearForm<TE>(LocalElementProductIntegrand<TE>()));
+          localizable_temporal_product.assemble();
+          current_data_["norm"][norm_id] = localizable_temporal_product.result();
+        } else
+          DUNE_THROW(XT::Common::Exceptions::wrong_input_given,
+                     "I do not know how to compute the temporal norm '" << temporal_norm_id << "'!");
+      }
+    }
+    // - estimates
+    std::vector<std::string> actual_estimates = self.filter(self.estimates(), only_these);
+    DUNE_THROW_IF(!actual_estimates.empty(),
+                  XT::Common::Exceptions::wrong_input_given,
+                  "I did not know how to compute the following estimates: " << actual_estimates);
+    return current_data_;
+  } // ... compute_on_current_refinement(...)
+
+protected:
+  static std::vector<double> time_points_from_vector_array(const XT::LA::ListVectorArray<V>& va)
+  {
+    std::vector<double> time_points;
+    time_points.reserve(va.length());
+    for (const auto& note : va.notes())
+      time_points.emplace_back(note.get("_t").at(0));
+    return time_points;
+  }
+
+  const double T_end_;
+  const size_t num_refinements_;
+  const size_t num_additional_refinements_for_reference_;
+  const std::function<void(const DiscreteBochnerFunction<V, GV, m>&, const std::string&)> visualize_;
+  size_t current_refinement_;
+  std::map<std::string, std::map<std::string, double>> current_data_;
+  std::unique_ptr<GP> current_grid_;
+  std::unique_ptr<S> current_space_;
+  std::unique_ptr<GP> reference_grid_;
+  std::unique_ptr<S> reference_space_;
+  std::unique_ptr<XT::LA::ListVectorArray<V>> current_solution_on_reference_grid_;
+  std::unique_ptr<XT::LA::ListVectorArray<V>> reference_solution_on_reference_grid_;
+}; // struct InstationaryEocStudy
+
+
+template <class V, class GV, size_t m>
+XT::LA::ListVectorArray<V> solve_instationary_system_explicit_euler(const DiscreteFunction<V, GV, m>& initial_values,
+                                                                    const OperatorInterface<V, GV, m>& spatial_op,
+                                                                    const double T_end,
+                                                                    const double dt)
+{
+  // initial values
+  XT::LA::ListVectorArray<V> solution(
+      spatial_op.source_space().mapper().size(), /*length=*/0, /*reserve=*/std::ceil(T_end / (dt)));
+  solution.append(initial_values.dofs().vector(), {"_t", 0.});
+  // timestepping
+  double time = 0.;
+  while (time < T_end + dt) {
+    time += dt;
+    const auto& u_n = solution.back().vector();
+    auto u_n_plus_one = u_n - spatial_op.apply(u_n, {{"_t", {time}}, {"_dt", {dt}}}) * dt;
+    solution.append(std::move(u_n_plus_one), {"_t", time});
+  }
+  return solution;
+} // ... solve_instationary_system_explicit_euler(...)
+
+
+} // namespace Tests
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_TEST_INSTATIONARY_EOCSTUDIES_BASE_HH

dune/gdt/test/instationary-eocstudy.hh

-// This file is part of the dune-gdt project:
-//   https://github.com/dune-community/dune-gdt
-// Copyright 2010-2018 dune-gdt developers and contributors. All rights reserved.
-// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
-//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
-//          with "runtime exception" (http://www.dune-project.org/license.html)
-// Authors:
-//   Felix Schindler (2016 - 2017)
-//   Rene Milk       (2016 - 2018)
-//   Tobias Leibner  (2016 - 2017)
-
-#ifndef DUNE_GDT_TEST_INSTATIONARY_EOCSTUDY_HH
-#define DUNE_GDT_TEST_INSTATIONARY_EOCSTUDY_HH
-
-#include <dune/xt/common/convergence-study.hh>
-#include <dune/xt/common/exceptions.hh>
-
-#include <dune/xt/grid/information.hh>
-
-#include <dune/gdt/projections.hh>
-#include <dune/gdt/prolongations.hh>
-
-namespace Dune {
-namespace GDT {
-namespace Test {
-
-
-template <class TestCaseImp, class DiscretizerImp>
-class InstationaryEocStudy : public XT::Common::ConvergenceStudy
-{
-  typedef XT::Common::ConvergenceStudy BaseType;
-
-protected:
-  typedef TestCaseImp TestCaseType;
-  typedef DiscretizerImp Discretizer;
-  typedef typename Discretizer::DiscretizationType DiscretizationType;
-  typedef typename TestCaseType::GridType GridType;
-  typedef typename Discretizer::ProblemType ProblemType;
-  typedef typename DiscretizationType::DiscreteSolutionType DiscreteSolutionType;
-  typedef typename DiscretizationType::DiscreteFunctionType DiscreteFunctionType;
-  typedef typename DiscreteSolutionType::key_type TimeFieldType;
-
-  typedef typename TestCaseType::InitialValueType InitialValueType;
-  typedef typename TestCaseType::LevelGridViewType GridLayerType;
-
-public:
-  InstationaryEocStudy(TestCaseType& test_case,
-                       const std::vector<std::string> only_these_norms = {},
-                       const std::string visualize_prefix = "transport_test")
-    : BaseType(only_these_norms)
-    , test_case_(test_case)
-    , current_refinement_(0)
-    , last_computed_refinement_(std::numeric_limits<size_t>::max())
-    , grid_widths_(num_refinements() + 1, -1.0)
-    , time_to_solution_(0)
-    , reference_solution_computed_(false)
-    , discrete_exact_solution_computed_(false)
-    , current_discretization_(nullptr)
-    , current_solution_on_level_(nullptr)
-    , reference_discretization_(nullptr)
-    , reference_solution_(nullptr)
-    , current_solution_(nullptr)
-    , discrete_exact_solution_(nullptr)
-    , visualize_prefix_(visualize_prefix)
-  {
-  }
-
-  virtual ~InstationaryEocStudy() = default;
-
-  virtual size_t num_refinements() const override final
-  {
-    return test_case_.num_refinements();
-  }
-
-  virtual std::vector<std::string> provided_norms() const override final
-  {
-    std::vector<std::string> ret = available_norms();
-    for (auto estimator : available_estimators()) {
-      if (is_norm(estimator))
-        DUNE_THROW(XT::Common::Exceptions::internal_error,
-                   "We do not want to handle the case that norms and estimators have the same name!");
-      ret.push_back(estimator);
-    }
-    return ret;
-  } // ... provided_norms(...)
-
-  virtual double norm_reference_solution(const std::string type) override final
-  {
-    if (is_norm(type)) {
-      compute_reference_solution();
-      assert(reference_solution_);
-      if (test_case_.provides_exact_solution()) {
-        compute_discrete_exact_solution_vector();
-        assert(discrete_exact_solution_);
-        return compute_norm(*discrete_exact_solution_, type);
-      } else {
-        return compute_norm(*reference_solution_, type);
-      }
-    } else
-      return 1.0;
-  } // ... norm_reference_solution(...)
-
-  virtual size_t current_num_DoFs() override final
-  {
-    assert(current_refinement_ <= num_refinements());
-    const int level = test_case_.level_of(current_refinement_);
-    return test_case_.grid().size(level, 0);
-  } // ... current_num_DoFs(...)
-
-  virtual size_t current_grid_size() const override final
-  {
-    assert(current_refinement_ <= num_refinements());
-    const int level = test_case_.level_of(current_refinement_);
-    return test_case_.grid().size(level, 0);
-  } // ... current_grid_size(...)
-
-  virtual double current_grid_width() override final
-  {
-    assert(current_refinement_ <= num_refinements());
-    if (grid_widths_[current_refinement_] < 0.0) {
-      const int level = test_case_.level_of(current_refinement_);
-      const auto grid_layer = test_case_.template layer<XT::Grid::Layers::level, XT::Grid::Backends::view>(level);
-      grid_widths_[current_refinement_] = XT::Grid::dimensions(grid_layer).entity_width.max();
-      assert(grid_widths_[current_refinement_] > 0.0);
-    }
-    return grid_widths_[current_refinement_];
-  } // ... current_grid_width(...)
-
-  virtual double compute_on_current_refinement() override final
-  {
-    if (current_refinement_ != last_computed_refinement_) {
-      assert(current_refinement_ <= num_refinements());
-      // compute solution
-      Timer timer;
-      current_discretization_ = XT::Common::make_unique<DiscretizationType>(Discretizer::discretize(
-          test_case_, test_case_, test_case_.level_of(current_refinement_), test_case_.periodic_directions()));
-      current_solution_on_level_ = XT::Common::make_unique<DiscreteSolutionType>(current_discretization_->solve());
-      time_to_solution_ = timer.elapsed();
-      // prolong to reference grid part
-      compute_reference_solution();
-      assert(reference_solution_);
-      if (!current_solution_)
-        current_solution_ = XT::Common::make_unique<DiscreteSolutionType>(*reference_solution_);
-      // time prolongation
-      DiscreteSolutionType time_prolongated_current_solution_on_level;
-      const auto time_prolongated_current_solution_on_level_it_end = time_prolongated_current_solution_on_level.end();
-      auto current_solution_on_level_it = current_solution_on_level_->begin();
-      const auto current_solution_on_level_it_last = --current_solution_on_level_->end();
-      const auto current_solution_it_end = current_solution_->end();
-      auto last_time = current_solution_->begin()->first;
-      for (auto current_solution_it = current_solution_->begin(); current_solution_it != current_solution_it_end;
-           ++current_solution_it) {
-        const auto time = current_solution_it->first;
-        const auto time_on_level = current_solution_on_level_it->first;
-        const auto inserted_it = time_prolongated_current_solution_on_level.emplace_hint(
-            time_prolongated_current_solution_on_level_it_end, time, current_solution_on_level_it->second);
-        if (time_on_level < time && current_solution_on_level_it != current_solution_on_level_it_last) {
-          // compute weighted average of the two values of current_solution_on_level_
-          inserted_it->second.vector() = (current_solution_on_level_it->second.vector() * (time_on_level - last_time)
-                                          + (++current_solution_on_level_it)->second.vector() * (time - time_on_level))
-                                         * (1.0 / (time - last_time));
-        }
-        last_time = time;
-      }
-      // spatial prolongation
-      for (auto current_solution_it = current_solution_->begin(); current_solution_it != current_solution_it_end;
-           ++current_solution_it)
-        prolong(time_prolongated_current_solution_on_level.at(current_solution_it->first),
-                (*current_solution_it).second);
-      last_computed_refinement_ = current_refinement_;
-      // visualize
-      if (!visualize_prefix_.empty()) {
-        size_t counter = 0;
-        for (auto current_solution_it = current_solution_->begin(); current_solution_it != current_solution_it_end;
-             ++current_solution_it, ++counter)
-          current_solution_it->second.visualize(visualize_prefix_ + "_solution_"
-                                                    + Dune::XT::Common::to_string(current_refinement_),
-                                                Dune::XT::Common::to_string(counter));
-      }
-    }
-    return time_to_solution_;
-  } // ... compute_on_current_refinement(...)
-
-  virtual double current_error_norm(const std::string type) override final
-  {
-    // get current solution
-    assert(current_refinement_ <= num_refinements());
-    compute_on_current_refinement();
-    assert(last_computed_refinement_ == current_refinement_);
-    if (is_norm(type)) {
-      assert(current_solution_);
-      compute_reference_solution();
-      assert(reference_discretization_);
-      // compute error
-      std::unique_ptr<DiscreteSolutionType> difference;
-      if (test_case_.provides_exact_solution()) {
-        compute_discrete_exact_solution_vector();
-        assert(discrete_exact_solution_);
-        difference = Dune::XT::Common::make_unique<DiscreteSolutionType>(*discrete_exact_solution_);
-      } else {
-        assert(reference_solution_);
-        difference = Dune::XT::Common::make_unique<DiscreteSolutionType>(*reference_solution_);
-      }
-      for (auto& difference_at_time : *difference) {
-        auto time = difference_at_time.first;
-        assert(current_solution_->count(time) && "Time steps must be the same");
-        difference_at_time.second.vector() = difference_at_time.second.vector() - current_solution_->at(time).vector();
-      }
-      return compute_norm(*difference, type);
-    } else {
-      assert(current_solution_on_level_);
-      return estimate(*current_solution_on_level_, type);
-    }
-  } // ... current_error_norm(...)
-
-  virtual void refine() override final
-  {
-    if (current_refinement_ <= num_refinements())
-      ++current_refinement_;
-  }
-
-  std::map<std::string, std::vector<double>> run(std::ostream& out, const bool print_timings = false)
-  {
-    return BaseType::run(false, out, print_timings);
-  }
-
-protected:
-  void compute_reference_solution()
-  {
-    if (!reference_solution_computed_) {
-      reference_discretization_ = XT::Common::make_unique<DiscretizationType>(Discretizer::discretize(
-          test_case_, test_case_, test_case_.reference_level(), test_case_.periodic_directions()));
-      reference_solution_ = XT::Common::make_unique<DiscreteSolutionType>(reference_discretization_->solve());
-      reference_solution_computed_ = true;
-      if (!visualize_prefix_.empty()) {
-        size_t counter = 0;
-        const auto reference_solution_it_end = reference_solution_->end();
-        for (auto reference_solution_it = reference_solution_->begin();
-             reference_solution_it != reference_solution_it_end;
-             ++reference_solution_it, ++counter)
-          reference_solution_it->second.visualize(visualize_prefix_ + "_reference"
-                                                      + Dune::XT::Common::to_string(current_refinement_),
-                                                  Dune::XT::Common::to_string(counter));
-      }
-    }
-  } // ... compute_reference_solution()
-
-  void compute_discrete_exact_solution_vector()
-  {
-    if (!discrete_exact_solution_computed_) {
-      discrete_exact_solution_ = Dune::XT::Common::make_unique<DiscreteSolutionType>();
-      compute_reference_solution();
-      const auto exact_solution = test_case_.exact_solution();
-      const auto reference_solution_it_end = reference_solution_->end();
-      for (auto reference_solution_it = reference_solution_->begin();
-           reference_solution_it != reference_solution_it_end;
-           ++reference_solution_it) {
-        const double time = reference_solution_it->first;
-        const auto inserted_it = discrete_exact_solution_->emplace_hint(
-            discrete_exact_solution_->end(), time, reference_solution_it->second);
-        project_l2(*exact_solution, inserted_it->second, 0, {"t", time});
-      }
-      if (!visualize_prefix_.empty()) {
-        size_t counter = 0;
-        const auto discrete_exact_solution_it_end = discrete_exact_solution_->end();
-        for (auto discrete_exact_solution_it = discrete_exact_solution_->begin();
-             discrete_exact_solution_it != discrete_exact_solution_it_end;
-             ++discrete_exact_solution_it, ++counter)
-          discrete_exact_solution_it->second.visualize(visualize_prefix_ + "_exact_solution"
-                                                           + Dune::XT::Common::to_string(current_refinement_),
-                                                       Dune::XT::Common::to_string(counter));
-      }
-      discrete_exact_solution_computed_ = true;
-    }
-  }
-
-  bool is_norm(const std::string type) const
-  {
-    const auto norms = available_norms();
-    return std::find(norms.begin(), norms.end(), type) != norms.end();
-  }
-
-  virtual std::vector<std::string> available_norms() const = 0;
-
-  virtual std::vector<std::string> available_estimators() const = 0;
-
-  virtual double estimate(const DiscreteSolutionType& solution, const std::string type) = 0;
-
-  virtual double compute_norm(const DiscreteSolutionType& solution, const std::string type) = 0;
-
-  TestCaseType& test_case_;
-  size_t current_refinement_;
-  size_t last_computed_refinement_;
-  mutable std::vector<double> grid_widths_;
-  double time_to_solution_;
-  bool reference_solution_computed_;
-  bool discrete_exact_solution_computed_;
-  std::unique_ptr<DiscretizationType> current_discretization_;
-  std::unique_ptr<DiscreteSolutionType> current_solution_on_level_;
-  std::unique_ptr<DiscretizationType> reference_discretization_;
-  std::unique_ptr<DiscreteSolutionType> reference_solution_;
-  std::unique_ptr<DiscreteSolutionType> current_solution_;
-  std::unique_ptr<DiscreteSolutionType> discrete_exact_solution_;
-  const std::string visualize_prefix_;
-}; // class InstationaryEocStudy
-
-
-} // namespace Tests
-} // namespace GDT
-} // namespace Dune
-
-#endif // DUNE_GDT_TEST_INSTATIONARY_EOCSTUDY_HH

dune/gdt/test/instationary-testcase.hh

-// This file is part of the dune-gdt project:
-//   https://github.com/dune-community/dune-gdt
-// Copyright 2010-2018 dune-gdt developers and contributors. All rights reserved.
-// License: Dual licensed as BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
-//      or  GPL-2.0+ (http://opensource.org/licenses/gpl-license)
-//          with "runtime exception" (http://www.dune-project.org/license.html)
-// Authors:
-//   Felix Schindler (2016 - 2017)
-//   Rene Milk       (2016 - 2018)
-//   Tobias Leibner  (2017)
-
-#ifndef DUNE_GDT_TEST_INSTATIONARY_TESTCASE_HH
-#define DUNE_GDT_TEST_INSTATIONARY_TESTCASE_HH
-
-#include <dune/xt/common/exceptions.hh>
-
-#include <dune/xt/grid/type_traits.hh>
-#include <dune/xt/grid/information.hh>
-#include <dune/xt/grid/gridprovider/eoc.hh>
-
-namespace Dune {
-namespace GDT {
-namespace Test {
-
-
-/**
- * \tparam ProblemType has to provide a type SolutionType which
- *         defines the type of the solution of the problem.
- * TODO: choose suitable SolutionType for Problems (provide Interface?)
- */
-template <class GridImp, class ProblemImp>
-class InstationaryTestCase : public XT::Grid::EOCGridProvider<GridImp>
-{
-  typedef XT::Grid::EOCGridProvider<GridImp> EocBaseType;
-
-public:
-  typedef ProblemImp ProblemType;
-  typedef typename ProblemType::InitialValueType InitialValueType;
-  typedef typename ProblemType::SolutionType SolutionType;
-
-public:
-  template <class... Args>
-  InstationaryTestCase(const double divide_t_end_by_this, Args&&... args)
-    : EocBaseType(std::forward<Args>(args)...)
-    , divide_t_end_by_this_(divide_t_end_by_this)
-    , zero_()
-  {
-  }
-
-  virtual ~InstationaryTestCase() = default;
-
-  virtual const ProblemType& problem() const = 0;
-
-  virtual void print_header(std::ostream& out = std::cout) const
-  {
-    out << "+===============================================================+\n"
-        << "|+=============================================================+|\n"
-        << "||  This is a GDT::Tests:InstationaryTestCase, please provide ||\n"
-        << "||  a meaningful message by implementing `print_header()`      ||\n"
-        << "|+=============================================================+|\n"
-        << "+===============================================================+" << std::endl;
-  }
-
-  virtual bool provides_exact_solution() const
-  {
-    return false;
-  }
-
-  virtual std::bitset<GridImp::dimension> periodic_directions() const
-  {
-    return std::bitset<GridImp::dimension>();
-  }
-
-  virtual double t_end() const
-  {
-    return problem().t_end() / divide_t_end_by_this_;
-  }
-
-  virtual const std::shared_ptr<const SolutionType> exact_solution() const
-  {
-    if (provides_exact_solution())
-      DUNE_THROW(XT::Common::Exceptions::you_have_to_implement_this,
-                 "If provides_exact_solution() is true, exact_solution() has to be implemented!");
-    else
-      DUNE_THROW(XT::Common::Exceptions::you_are_using_this_wrong,
-                 "Do not call exact_solution() if provides_exact_solution() is false!");
-    return zero_;
-  }
-
-private:
-  const double divide_t_end_by_this_;
-  const std::shared_ptr<const SolutionType> zero_;
-}; // class InstationaryTestCase
-
-
-} // namespace Tests
-} // namespace GDT
-} // namespace Dune
-
-#endif // DUNE_GDT_TEST_INSTATIONARY_TESTCASE_HH