From 401c98cf1b376c8821634d15678a4feef868ed64 Mon Sep 17 00:00:00 2001
From: Felix Schindler <felix.schindler@wwu.de>
Date: Sun, 5 Aug 2018 17:56:38 +0200
Subject: [PATCH] [test] complete rewrite and modernize of InstationaryEocStudy
---
dune/gdt/test/instationary-eocstudies/base.hh | 338 ++++++++++++++++++
dune/gdt/test/instationary-eocstudy.hh | 312 ----------------
dune/gdt/test/instationary-testcase.hh | 100 ------
3 files changed, 338 insertions(+), 412 deletions(-)
create mode 100644 dune/gdt/test/instationary-eocstudies/base.hh
delete mode 100644 dune/gdt/test/instationary-eocstudy.hh
delete mode 100644 dune/gdt/test/instationary-testcase.hh
diff --git a/dune/gdt/test/instationary-eocstudies/base.hh b/dune/gdt/test/instationary-eocstudies/base.hh
new file mode 100644
index 000000000..372d5967f
--- /dev/null
+++ b/dune/gdt/test/instationary-eocstudies/base.hh
@@ -0,0 +1,338 @@
+// 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
diff --git a/dune/gdt/test/instationary-eocstudy.hh b/dune/gdt/test/instationary-eocstudy.hh
deleted file mode 100644
index dcfeb37cb..000000000
--- a/dune/gdt/test/instationary-eocstudy.hh
+++ /dev/null
@@ -1,312 +0,0 @@
-// 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
diff --git a/dune/gdt/test/instationary-testcase.hh b/dune/gdt/test/instationary-testcase.hh
deleted file mode 100644
index c7a2bc8b3..000000000
--- a/dune/gdt/test/instationary-testcase.hh
+++ /dev/null
@@ -1,100 +0,0 @@
-// 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
--
GitLab