diff --git a/dune/gdt/local/numerical-fluxes/kinetic.hh b/dune/gdt/local/numerical-fluxes/kinetic.hh
index c0c7962642ac314d354f9fc3dbd803fd0e324f34..9c1f339755ac231a3973934538b0518b205ec58d 100644
--- a/dune/gdt/local/numerical-fluxes/kinetic.hh
+++ b/dune/gdt/local/numerical-fluxes/kinetic.hh
@@ -22,7 +22,7 @@ namespace Dune {
namespace GDT {
-template <class GV, class MomentBasis>
+template <class GV, class MomentBasis, class EntropyFluxImp = EntropyBasedFluxFunction<GV, MomentBasis>>
class NumericalKineticFlux
: public NumericalFluxInterface<XT::Grid::extract_intersection_t<GV>,
MomentBasis::dimFlux,
@@ -43,7 +43,7 @@ class NumericalKineticFlux
using I = XT::Grid::extract_intersection_t<GV>;
using ThisType = NumericalKineticFlux;
using BaseType = NumericalFluxInterface<I, d, m, R>;
- using EntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
+ using EntropyFluxType = EntropyFluxImp;
using SparseMatrixType = typename XT::LA::CommonSparseMatrix<R>;
public:
@@ -83,13 +83,22 @@ public:
// find direction of unit outer normal (we assume an axis-aligned cube grid)
size_t direction = intersection().indexInInside() / 2;
if (dynamic_cast<const EntropyFluxType*>(&flux()) != nullptr) {
- return dynamic_cast<const EntropyFluxType*>(&flux())->evaluate_kinetic_flux(
+ auto ret = dynamic_cast<const EntropyFluxType*>(&flux())->evaluate_kinetic_flux(
intersection().inside(),
intersection().neighbor() ? intersection().outside() : intersection().inside(),
u,
v,
n,
direction);
+ for (auto&& entry : ret)
+ if (std::isnan(entry) || std::isinf(entry)) {
+ // std::cout << XT::Common::to_string(ret) << std::endl;
+ // std::cout << XT::Common::to_string(u) << std::endl;
+ // std::cout << XT::Common::to_string(v) << std::endl;
+ // std::cout << this->intersection().geometry().center() << std::endl;
+ DUNE_THROW(Dune::MathError, "NaN or inf in kinetic flux");
+ }
+ return ret;
} else {
static const auto flux_matrices = initialize_flux_matrices(basis_);
StateType ret(0);
diff --git a/dune/gdt/momentmodels/entropyflux.hh b/dune/gdt/momentmodels/entropyflux.hh
index 02a5e54299be08362537565ad8d9c158ceb9f7b2..ca237ecb29b6c04a3fb0bf85db35d46ca27b2a4e 100644
--- a/dune/gdt/momentmodels/entropyflux.hh
+++ b/dune/gdt/momentmodels/entropyflux.hh
@@ -161,7 +161,8 @@ public:
, use_entity_cache_(true)
, entity_caches_(index_set_.size(0), LocalCacheType(cache_size))
, mutexes_(index_set_.size(0))
- , implementation_(basis_functions, tau, epsilon_gamma, chi, xi, r_sequence, k_0, k_max, epsilon)
+ , implementation_(std::make_shared<ImplementationType>(
+ basis_functions, tau, epsilon_gamma, chi, xi, r_sequence, k_0, k_max, epsilon))
{}
void enable_thread_cache()
@@ -303,13 +304,13 @@ public:
virtual std::unique_ptr<LocalFunctionType> local_function() const override final
{
return std::make_unique<Localfunction>(
- index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, implementation_);
+ index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, *implementation_);
}
virtual std::unique_ptr<Localfunction> derived_local_function() const
{
return std::make_unique<Localfunction>(
- index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, implementation_);
+ index_set_, entity_caches_, use_thread_cache_, use_entity_cache_, mutexes_, *implementation_);
}
StateType evaluate_kinetic_flux(const E& inside_entity,
@@ -325,21 +326,20 @@ public:
const auto alpha_i = local_func->get_alpha(u_i, true)->first;
local_func->bind(outside_entity);
const auto alpha_j = local_func->get_alpha(u_j, true)->first;
- return implementation_.evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
+ return implementation_->evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
} // StateType evaluate_kinetic_flux(...)
const MomentBasis& basis_functions() const
{
- return implementation_.basis_functions();
+ return implementation_->basis_functions();
}
-private:
const IndexSetType& index_set_;
bool use_thread_cache_;
bool use_entity_cache_;
mutable std::vector<LocalCacheType> entity_caches_;
mutable std::vector<std::mutex> mutexes_;
- ImplementationType implementation_;
+ std::shared_ptr<ImplementationType> implementation_;
};
template <class GridViewImp, class MomentBasisImp>
diff --git a/dune/gdt/momentmodels/entropyflux_implementations.hh b/dune/gdt/momentmodels/entropyflux_implementations.hh
index 2b8429f95fadec68675265deebfcd60685c07b43..5d3fdbc3052db4ea68092642b77154a05ff3baaa 100644
--- a/dune/gdt/momentmodels/entropyflux_implementations.hh
+++ b/dune/gdt/momentmodels/entropyflux_implementations.hh
@@ -105,6 +105,7 @@ public:
using DynamicRangeType = DynamicVector<RangeFieldType>;
using BasisValuesMatrixType = XT::LA::CommonDenseMatrix<RangeFieldType>;
using AlphaReturnType = std::pair<VectorType, std::pair<DomainType, RangeFieldType>>;
+ using QuadraturePointsValuesType = std::vector<RangeFieldType>;
explicit EntropyBasedFluxImplementationUnspecializedBase(const MomentBasis& basis_functions,
const RangeFieldType tau,
@@ -177,7 +178,7 @@ public:
return evaluate_with_alpha(alpha);
}
- virtual RangeReturnType evaluate_with_alpha(const VectorType& alpha) const
+ virtual RangeReturnType evaluate_with_alpha(const DomainType& alpha) const
{
RangeReturnType ret(0.);
auto& work_vecs = working_storage();
@@ -202,7 +203,7 @@ public:
jacobian_with_alpha(alpha, result);
}
- virtual void jacobian_with_alpha(const VectorType& alpha, DynamicDerivativeRangeType& result) const
+ virtual void jacobian_with_alpha(const DomainType& alpha, DynamicDerivativeRangeType& result) const
{
thread_local auto H = XT::Common::make_unique<MatrixType>();
calculate_hessian(alpha, M_, *H);
@@ -223,8 +224,8 @@ public:
evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
} // DomainType evaluate_kinetic_flux(...)
- DomainType evaluate_kinetic_flux_with_alphas(const VectorType& alpha_i,
- const VectorType& alpha_j,
+ DomainType evaluate_kinetic_flux_with_alphas(const DomainType& alpha_i,
+ const DomainType& alpha_j,
const FluxDomainType& n_ij,
const size_t dd) const
@@ -254,9 +255,8 @@ public:
// returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
virtual std::unique_ptr<AlphaReturnType>
- get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const = 0;
+ get_alpha(const DomainType& u, const DomainType& alpha_in, const bool regularize) const = 0;
-protected:
// get permutation instead of sorting directly to be able to sort two vectors the same way
// see
// https://stackoverflow.com/questions/17074324/how-can-i-sort-two-vectors-in-the-same-way-with-criteria-that-uses-only-one-of
@@ -439,7 +439,7 @@ protected:
return work_vec;
}
- void calculate_scalar_products(const VectorType& beta_in,
+ void calculate_scalar_products(const DomainType& beta_in,
const BasisValuesMatrixType& M,
std::vector<RangeFieldType>& scalar_products) const
{
@@ -472,7 +472,7 @@ protected:
}
// calculate ret = \int (exp(beta_in * m))
- RangeFieldType calculate_scalar_integral(const VectorType& beta_in, const BasisValuesMatrixType& M) const
+ RangeFieldType calculate_scalar_integral(const DomainType& beta_in, const BasisValuesMatrixType& M) const
{
auto& work_vec = working_storage();
calculate_scalar_products(beta_in, M, work_vec);
@@ -481,10 +481,10 @@ protected:
}
// calculate ret = \int (m1 exp(beta_in * m2))
- void calculate_vector_integral(const VectorType& beta_in,
+ void calculate_vector_integral(const DomainType& beta_in,
const BasisValuesMatrixType& M1,
const BasisValuesMatrixType& M2,
- VectorType& ret,
+ DomainType& ret,
bool same_beta = false,
bool only_first_component = false) const
{
@@ -571,7 +571,17 @@ protected:
} // ii
} // void calculate_A_Binv(...)
- void calculate_hessian(const VectorType& alpha,
+ void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ {
+ thread_local auto H = XT::Common::make_unique<MatrixType>();
+ calculate_hessian(alpha, M_, *H);
+ XT::LA::cholesky(*H);
+ thread_local DomainType tmp_vec;
+ XT::LA::solve_lower_triangular(*H, tmp_vec, u);
+ XT::LA::solve_lower_triangular_transposed(*H, Hinv_u, tmp_vec);
+ }
+
+ void calculate_hessian(const DomainType& alpha,
const BasisValuesMatrixType& M,
MatrixType& H,
const bool use_work_vec_data = false) const
@@ -629,12 +639,12 @@ protected:
J_dd.set_entry(mm, nn, J_dd.get_entry(nn, mm));
} // void calculate_J(...)
- void change_basis(const VectorType& beta_in,
- VectorType& v_k,
+ void change_basis(const DomainType& beta_in,
+ DomainType& v_k,
BasisValuesMatrixType& P_k,
MatrixType& T_k,
- VectorType& g_k,
- VectorType& beta_out,
+ DomainType& g_k,
+ DomainType& beta_out,
MatrixType& H) const
{
calculate_hessian(beta_in, P_k, H);
@@ -683,12 +693,15 @@ class EntropyBasedFluxImplementation : public EntropyBasedFluxImplementationUnsp
using ThisType = EntropyBasedFluxImplementation;
public:
+ using BaseType::basis_dimRange;
using BaseType::dimFlux;
using typename BaseType::AlphaReturnType;
+ using typename BaseType::BasisDomainType;
using typename BaseType::BasisValuesMatrixType;
using typename BaseType::DomainType;
using typename BaseType::MatrixType;
using typename BaseType::MomentBasis;
+ using typename BaseType::QuadraturePointsValuesType;
using typename BaseType::RangeFieldType;
using typename BaseType::VectorType;
@@ -707,9 +720,119 @@ public:
XT::LA::eye_matrix(*T_minus_one_);
}
+ 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
+ {
+ DomainType ret;
+ calculate_vector_integral(alpha, M_, M_, ret);
+ return ret;
+ }
+
+ static RangeFieldType minmod(const RangeFieldType& first_slope, const RangeFieldType& second_slope)
+ {
+ RangeFieldType ret = 0.;
+ if (std::signbit(first_slope) == std::signbit(second_slope)) // check for equal sign
+ ret = std::abs(first_slope) < std::abs(second_slope) ? first_slope : second_slope;
+ return ret;
+ }
+
+ static RangeFieldType maxmod(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ {
+ RangeFieldType ret(0.);
+ if (std::signbit(first_slope) == std::signbit(second_slope)) // check for equal sign
+ ret = (std::abs(first_slope) < std::abs(second_slope)) ? second_slope : first_slope;
+ return ret;
+ }
+
+ static RangeFieldType superbee(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ {
+ const RangeFieldType first_minmod = minmod(first_slope, second_slope * 2.);
+ const RangeFieldType second_minmod = minmod(first_slope * 2., second_slope);
+ return maxmod(first_minmod, second_minmod);
+ }
+
+ template <class IntersectionType>
+ DomainType calculate_boundary_flux(const DomainType& alpha, const IntersectionType& intersection)
+ {
+ // evaluate exp(alpha^T b(v_i)) at all quadratures points v_i
+ auto& work_vecs = this->working_storage();
+ this->calculate_scalar_products(alpha, M_, work_vecs);
+ this->apply_exponential(work_vecs);
+ const auto dd = intersection.indexInInside() / 2;
+ const auto& n_ij = intersection.centerUnitOuterNormal();
+ DomainType ret(0.);
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll) {
+ const auto position = quad_points_[ll][dd];
+ if (position * n_ij[dd] < 0.) {
+ RangeFieldType factor = work_vecs[ll] * quad_weights_[ll] * position;
+ const auto* basis_ll = &(M_.get_entry_ref(ll, 0.));
+ for (size_t ii = 0; ii < basis_dimRange; ++ii)
+ ret[ii] += basis_ll[ii] * factor;
+ }
+ } // ll
+ return ret;
+ }
+
+ // TODO: Calculates exp(alpha^T b) 2*dim times for each alpha (as it is contained in 2*dim stencils) at the moment.
+ template <class FluxesMapType>
+ void calculate_minmod_density_reconstruction(const FieldVector<DomainType, 3>& alphas,
+ FluxesMapType& flux_values,
+ const size_t dd) const
+ {
+ // evaluate exp(alpha^T b(v_i)) at all quadratures points v_i for all three alphas
+ thread_local XT::Common::FieldVector<QuadraturePointsValuesType, 3> density_evaluations(
+ QuadraturePointsValuesType(quad_points_.size()));
+ thread_local XT::Common::FieldVector<QuadraturePointsValuesType, 2> reconstructed_values(
+ QuadraturePointsValuesType(quad_points_.size()));
+ for (size_t ii = 0; ii < 3; ++ii) {
+ this->calculate_scalar_products(alphas[ii], M_, density_evaluations[ii]);
+ this->apply_exponential(density_evaluations[ii]);
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll) {
+ if (std::isnan(density_evaluations[ii][ll]) || std::isnan(density_evaluations[ii][ll])) {
+ std::cout << "Wrong val: " << density_evaluations[ii][ll] << " for alpha " << alphas[ii] << std::endl;
+ DUNE_THROW(Dune::MathError, "");
+ }
+ }
+ }
+ // get left and right reconstructed values for each quadrature point v_i
+ auto& vals_left = reconstructed_values[0];
+ auto& vals_right = reconstructed_values[1];
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll) {
+ const auto slope = minmod(density_evaluations[1][ll] - density_evaluations[0][ll],
+ density_evaluations[2][ll] - density_evaluations[1][ll]);
+ // const auto slope = superbee(density_evaluations[1][ll] - density_evaluations[0][ll],
+ // density_evaluations[2][ll] - density_evaluations[1][ll]);
+ vals_left[ll] = density_evaluations[1][ll] - 0.5 * slope;
+ vals_right[ll] = density_evaluations[1][ll] + 0.5 * slope;
+ } // ll
+
+ BasisDomainType coord(0.5);
+ coord[dd] = 0;
+ auto& left_flux_value = flux_values[coord];
+ coord[dd] = 1;
+ auto& right_flux_value = flux_values[coord];
+ right_flux_value = left_flux_value = DomainType(0.);
+
+ for (size_t ll = 0; ll < quad_points_.size(); ++ll) {
+ const auto position = quad_points_[ll][dd];
+ RangeFieldType factor = position > 0. ? vals_right[ll] : vals_left[ll];
+ factor *= quad_weights_[ll] * position;
+ auto& val = position > 0. ? right_flux_value : left_flux_value;
+ const auto* basis_ll = &(M_.get_entry_ref(ll, 0.));
+ for (size_t ii = 0; ii < basis_dimRange; ++ii)
+ val[ii] += basis_ll[ii] * factor;
+ } // ll
+ // std::cout << "right_flux_value, left_flux_value = " << right_flux_value << ", " << left_flux_value <<
+ // std::endl;
+ } // void calculate_minmod_density_reconstruction(...)
+
// returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
virtual std::unique_ptr<AlphaReturnType>
- get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const override final
+ get_alpha(const DomainType& u, const DomainType& alpha_in, const bool regularize) const override final
{
auto ret = std::make_unique<AlphaReturnType>();
@@ -841,12 +964,12 @@ private:
using BaseType::calculate_vector_integral;
using BaseType::get_isotropic_alpha;
- void change_basis(const VectorType& beta_in,
- VectorType& v_k,
+ void change_basis(const DomainType& beta_in,
+ DomainType& v_k,
BasisValuesMatrixType& P_k,
MatrixType& T_k,
- VectorType& g_k,
- VectorType& beta_out,
+ DomainType& g_k,
+ DomainType& beta_out,
MatrixType& H) const
{
calculate_hessian(beta_in, P_k, H);
@@ -913,6 +1036,16 @@ public:
: BaseType(basis_functions, tau, epsilon_gamma, chi, xi, r_sequence, k_0, k_max, epsilon)
{}
+ 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_vector_integral(alpha, M_, M_);
+ }
+
// returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
virtual std::unique_ptr<AlphaReturnType>
get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const override final
@@ -1233,6 +1366,18 @@ public:
return basis_functions_;
}
+ 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
+ {
+ BlockVectorType u_block;
+ calculate_vector_integral(alpha, M_, M_, u_block);
+ return u_block;
+ }
+
std::unique_ptr<AlphaReturnType>
get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const
{
@@ -1372,7 +1517,6 @@ public:
return ret;
}
-private:
// temporary vectors to store inner products and exponentials
TemporaryVectorsType& working_storage() const
{
@@ -1628,6 +1772,26 @@ private:
} // jj
} // void calculate_A_Binv(...)
+ void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ {
+ thread_local auto H = XT::Common::make_unique<BlockMatrixType>();
+ calculate_hessian(alpha, M_, *H);
+ for (size_t jj = 0; jj < num_blocks; ++jj)
+ XT::LA::cholesky(H->block(jj));
+ FieldVector<RangeFieldType, block_size> tmp_vec;
+ FieldVector<RangeFieldType, block_size> block_u;
+ for (size_t jj = 0; jj < num_blocks; ++jj) {
+ // calculate C = A (L^T)^{-1}
+ const auto offset = block_size * jj;
+ for (size_t ii = 0; ii < block_size; ++ii)
+ block_u[ii] = u[offset + ii];
+ XT::LA::solve_lower_triangular(H->block(jj), tmp_vec, block_u);
+ XT::LA::solve_lower_triangular_transposed(H->block(jj), block_u, tmp_vec);
+ for (size_t ii = 0; ii < block_size; ++ii)
+ Hinv_u[offset + ii] = block_u[ii];
+ } // jj
+ }
+
void calculate_hessian(const BlockVectorType& alpha, const BasisValuesMatrixType& M, BlockMatrixType& H) const
{
auto& work_vec = working_storage();
@@ -1858,6 +2022,13 @@ public:
return evaluate_with_alpha(alpha);
}
+ RangeReturnType evaluate_with_alpha(const DomainType& alpha) const
+ {
+ VectorType alpha_vec(basis_dimRange, 0., 0);
+ std::copy(alpha.begin(), alpha.end(), alpha_vec.begin());
+ return evaluate_with_alpha(alpha_vec);
+ }
+
virtual RangeReturnType evaluate_with_alpha(const VectorType& alpha) const
{
RangeReturnType ret(0.);
@@ -1893,7 +2064,14 @@ public:
jacobian_with_alpha(alpha, result);
}
- virtual void jacobian_with_alpha(const VectorType& alpha, DynamicDerivativeRangeType& result) const
+ void jacobian_with_alpha(const DomainType& alpha, DynamicDerivativeRangeType& result) const
+ {
+ VectorType alpha_vec(basis_dimRange, 0., 0);
+ std::copy(alpha.begin(), alpha.end(), alpha_vec.begin());
+ jacobian_with_alpha(alpha_vec, result);
+ }
+
+ void jacobian_with_alpha(const VectorType& alpha, DynamicDerivativeRangeType& result) const
{
thread_local SparseMatrixType H(basis_dimRange, basis_dimRange, pattern_, 0);
thread_local SparseMatrixType J(basis_dimRange, basis_dimRange, pattern_, 0);
@@ -1912,6 +2090,18 @@ public:
evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
} // DomainType evaluate_kinetic_flux(...)
+ DomainType evaluate_kinetic_flux_with_alphas(const DomainType& alpha_i,
+ const DomainType& alpha_j,
+ const FluxDomainType& n_ij,
+ const size_t dd) const
+ {
+ VectorType alpha_i_vec(basis_dimRange);
+ VectorType alpha_j_vec(basis_dimRange);
+ std::copy(alpha_i.begin(), alpha_i.end(), alpha_i_vec.begin());
+ std::copy(alpha_j.begin(), alpha_j.end(), alpha_j_vec.begin());
+ return evaluate_kinetic_flux_with_alphas(alpha_i_vec, alpha_j_vec, n_ij, dd);
+ }
+
DomainType evaluate_kinetic_flux_with_alphas(const VectorType& alpha_i,
const VectorType& alpha_j,
const FluxDomainType& n_ij,
@@ -1951,6 +2141,21 @@ public:
return basis_functions_;
}
+ 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
+ {
+ VectorType u_vec(basis_dimRange), alpha_vec(basis_dimRange);
+ std::copy(alpha.begin(), alpha.end(), alpha_vec.begin());
+ calculate_u(alpha_vec, u_vec);
+ DomainType u;
+ std::copy(u_vec.begin(), u_vec.end(), u.begin());
+ return u;
+ }
+
std::unique_ptr<AlphaReturnType>
get_alpha(const DomainType& u, const VectorType& alpha_in, const bool regularize) const
{
@@ -2087,7 +2292,6 @@ public:
return ret;
} // ... get_alpha(...)
-private:
static bool is_realizable(const VectorType& u)
{
for (const auto& u_i : u)
@@ -2108,7 +2312,6 @@ private:
return work_vecs;
}
-private:
// calculates ret = J H^{-1}. H is assumed to be symmetric positive definite, which gives ret^T = H^{-T} J^T =
// H^{-1} J^T, so we just have to solve y = H^{-1} x for each row x of J
void calculate_J_Hinv(SparseMatrixType& J, const SparseMatrixType& H, DynamicRowDerivativeRangeType& ret) const
@@ -2140,6 +2343,29 @@ private:
} // ii
} // void calculate_J_Hinv(...)
+ void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ {
+ thread_local SparseMatrixType H(basis_dimRange, basis_dimRange, pattern_, 0);
+ thread_local VectorType alpha_vec(basis_dimRange, 0., 0), u_vec(basis_dimRange, 0., 0),
+ Hinv_u_vec(basis_dimRange, 0., 0);
+ std::copy(alpha.begin(), alpha.end(), alpha_vec.begin());
+ std::copy(u.begin(), u.end(), u_vec.begin());
+ calculate_hessian(alpha_vec, M_, H);
+# if HAVE_EIGEN
+ typedef ::Eigen::SparseMatrix<RangeFieldType, ::Eigen::ColMajor> ColMajorBackendType;
+ ColMajorBackendType colmajor_copy(H.backend());
+ typedef ::Eigen::SimplicialLDLT<ColMajorBackendType> SolverType;
+ SolverType solver;
+ solver.analyzePattern(colmajor_copy);
+ solver.factorize(colmajor_copy);
+ Hinv_u_vec.backend() = solver.solve(u_vec.backend());
+ std::copy(Hinv_u_vec.begin(), Hinv_u_vec.end(), Hinv_u.begin());
+# else // HAVE_EIGEN
+ auto solver = XT::LA::make_solver(H);
+ solver.apply(u, Hinv_u);
+# endif
+ } // void apply_inverse_hessian(..)
+
RangeFieldType calculate_f(const VectorType& alpha, const VectorType& v) const
{
RangeFieldType ret(0.);
@@ -2598,6 +2824,15 @@ public:
return ret;
} // DomainType evaluate_kinetic_flux(...)
+ 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_vector_integral(alpha, M_, M_);
+ }
// returns (alpha, (actual_u, r)), where r is the regularization parameter and actual_u the regularized u
std::unique_ptr<AlphaReturnType>
@@ -3167,6 +3402,120 @@ public:
return ret;
} // ... evaluate_kinetic_flux(...)
+ 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
+ {
+ DomainType ret;
+ calculate_u(alpha, ret);
+ return ret;
+ }
+
+ static RangeFieldType minmod(const RangeFieldType& first_slope, const RangeFieldType& second_slope)
+ {
+ RangeFieldType ret = 0.;
+ if (std::signbit(first_slope) == std::signbit(second_slope)) // check for equal sign
+ ret = std::abs(first_slope) < std::abs(second_slope) ? first_slope : second_slope;
+ return ret;
+ }
+
+ static RangeFieldType maxmod(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ {
+ RangeFieldType ret(0.);
+ if (std::signbit(first_slope) == std::signbit(second_slope)) // check for equal sign
+ ret = (std::abs(first_slope) < std::abs(second_slope)) ? second_slope : first_slope;
+ return ret;
+ }
+
+ static RangeFieldType superbee(const RangeFieldType first_slope, const RangeFieldType second_slope)
+ {
+ const RangeFieldType first_minmod = minmod(first_slope, second_slope * 2.);
+ const RangeFieldType second_minmod = minmod(first_slope * 2., second_slope);
+ return maxmod(first_minmod, second_minmod);
+ }
+
+ template <class IntersectionType>
+ DomainType calculate_boundary_flux(const DomainType& alpha, const IntersectionType& intersection)
+ {
+ // evaluate exp(alpha^T b(v_i)) at all quadratures points v_i
+ const auto dd = intersection.indexInInside() / 2;
+ const auto& n_ij = intersection.centerUnitOuterNormal();
+ DomainType ret(0.);
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ for (size_t ll = 0; ll < quad_points_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll];
+ if (position * n_ij[dd] < 0.) {
+ RangeFieldType factor =
+ std::exp(alpha[jj] * M_[jj][ll][0] + alpha[jj + 1] * M_[jj][ll][1]) * quad_weights_[jj][ll] * position;
+ const auto& basis_ll = M_[jj][ll];
+ for (size_t mm = 0; mm < 2; ++mm)
+ ret[jj + mm] += basis_ll[mm] * factor;
+ }
+ } // ll
+ } // jj
+ return ret;
+ }
+
+ // TODO: Calculates exp(alpha^T b) 2*dim times for each alpha (as it is contained in 2*dim stencils) at the moment.
+ template <class FluxesMapType>
+ void calculate_minmod_density_reconstruction(const FieldVector<DomainType, 3>& alphas,
+ FluxesMapType& flux_values,
+ const size_t dd) const
+ {
+ // evaluate exp(alpha^T b(v_i)) at all quadratures points v_i for all three alphas
+ thread_local XT::Common::FieldVector<std::vector<RangeFieldType>, 3> density_evaluations(
+ std::vector<RangeFieldType>(quad_points_[0].size()));
+ thread_local XT::Common::FieldVector<std::vector<RangeFieldType>, 2> reconstructed_values(
+ std::vector<RangeFieldType>(quad_points_[0].size()));
+ thread_local XT::Common::FieldVector<LocalVectorType, 3> local_alphas;
+
+ // get flux storage
+ BasisDomainType coord(0.5);
+ coord[dd] = 0;
+ auto& left_flux_value = flux_values[coord];
+ coord[dd] = 1;
+ auto& right_flux_value = flux_values[coord];
+ right_flux_value = left_flux_value = DomainType(0.);
+
+ for (size_t jj = 0; jj < num_intervals; ++jj) {
+ // get local alphas
+ for (size_t ii = 0; ii < 3; ++ii)
+ for (size_t mm = 0; mm < 2; ++mm)
+ local_alphas[ii][mm] = alphas[ii][jj + mm];
+
+ // evaluate densitys
+ for (size_t ii = 0; ii < 3; ++ii)
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll)
+ density_evaluations[ii][ll] = std::exp(local_alphas[ii] * M_[jj][ll]);
+
+ // reconstruct densities
+ auto& vals_left = reconstructed_values[0];
+ auto& vals_right = reconstructed_values[1];
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto slope = minmod(density_evaluations[1][ll] - density_evaluations[0][ll],
+ density_evaluations[2][ll] - density_evaluations[1][ll]);
+ // const auto slope = superbee(density_evaluations[1][ll] - density_evaluations[0][ll],
+ // density_evaluations[2][ll] - density_evaluations[1][ll]);
+ vals_left[ll] = density_evaluations[1][ll] - 0.5 * slope;
+ vals_right[ll] = density_evaluations[1][ll] + 0.5 * slope;
+ } // ll (quad points)
+
+ // calculate fluxes
+ for (size_t ll = 0; ll < quad_weights_[jj].size(); ++ll) {
+ const auto position = quad_points_[jj][ll];
+ RangeFieldType factor = position > 0. ? vals_right[ll] : vals_left[ll];
+ factor *= quad_weights_[jj][ll] * position;
+ auto& val = position > 0. ? right_flux_value : left_flux_value;
+ const auto& basis_ll = M_[jj][ll];
+ for (size_t mm = 0; mm < 2; ++mm)
+ val[jj + mm] += basis_ll[mm] * factor;
+ } // ll (quad points)
+ } // jj (intervals)
+ } // void calculate_minmod_density_reconstruction(...)
+
// 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
@@ -3277,7 +3626,6 @@ public:
return basis_functions_;
}
-private:
static bool is_realizable(const DomainType& u)
{
for (const auto& u_i : u)
@@ -3309,7 +3657,7 @@ private:
return ret;
} // void calculate_u(...)
- void calculate_u(const VectorType& alpha, VectorType& u) const
+ void calculate_u(const DomainType& alpha, DomainType& u) const
{
std::fill(u.begin(), u.end(), 0.);
LocalVectorType local_alpha;
@@ -3331,7 +3679,7 @@ private:
g_k -= v;
}
- void calculate_hessian(const VectorType& alpha,
+ void calculate_hessian(const DomainType& alpha,
const BasisValuesMatrixType& M,
VectorType& H_diag,
FieldVector<RangeFieldType, dimRange - 1>& H_subdiag) const
@@ -3377,6 +3725,21 @@ private:
} // jj (intervals)
} // void calculate_J(...)
+ void apply_inverse_hessian(const DomainType& alpha, const DomainType& u, DomainType& Hinv_u) const
+ {
+ thread_local VectorType H_diag;
+ thread_local FieldVector<RangeFieldType, dimRange - 1> H_subdiag;
+ calculate_hessian(alpha, M_, H_diag, H_subdiag);
+ // std::cout << "H_diag: " << H_diag << ", H_subdiag: " << H_subdiag << std::endl;
+ // factorize H = LDL^T, where L is unit lower bidiagonal and D is diagonal
+ // H_diag is overwritten by the diagonal elements of D
+ // H_subdiag is overwritten by the subdiagonal elements of L
+ XT::LA::tridiagonal_ldlt(H_diag, H_subdiag);
+ // Solve H ret = u
+ Hinv_u = u;
+ XT::LA::solve_tridiagonal_ldlt_factorized(H_diag, H_subdiag, Hinv_u);
+ }
+
// calculates ret = J H^{-1}. Both J and H are symmetric tridiagonal, H is positive definite.
static void calculate_J_Hinv(DynamicRowDerivativeRangeType& ret,
const VectorType& J_diag,
diff --git a/dune/gdt/momentmodels/entropyflux_kineticcoords.hh b/dune/gdt/momentmodels/entropyflux_kineticcoords.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4dafa99a3ef22cd762e1e9904bf4238eac80586e
--- /dev/null
+++ b/dune/gdt/momentmodels/entropyflux_kineticcoords.hh
@@ -0,0 +1,195 @@
+// 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:
+// Tobias Leibner (2019)
+
+#ifndef DUNE_GDT_LOCAL_FLUXES_ENTROPYBASED_KINETICCOORDS_HH
+#define DUNE_GDT_LOCAL_FLUXES_ENTROPYBASED_KINETICCOORDS_HH
+
+#include <list>
+#include <memory>
+
+#include <dune/xt/common/float_cmp.hh>
+#include <dune/xt/common/vector_less.hh>
+
+#include <dune/xt/functions/interfaces/flux-function.hh>
+
+#include <dune/gdt/momentmodels/basisfunctions.hh>
+#include <dune/gdt/momentmodels/entropyflux_implementations.hh>
+#include <dune/gdt/momentmodels/entropyflux.hh>
+
+namespace Dune {
+namespace GDT {
+
+
+template <class GridViewImp, class MomentBasisImp>
+class EntropyBasedFluxEntropyCoordsFunction
+ : public XT::Functions::FluxFunctionInterface<XT::Grid::extract_entity_t<GridViewImp>,
+ MomentBasisImp::dimRange,
+ MomentBasisImp::dimFlux,
+ MomentBasisImp::dimRange,
+ typename MomentBasisImp::R>
+{
+ using BaseType = typename XT::Functions::FluxFunctionInterface<XT::Grid::extract_entity_t<GridViewImp>,
+ MomentBasisImp::dimRange,
+ MomentBasisImp::dimFlux,
+ MomentBasisImp::dimRange,
+ typename MomentBasisImp::R>;
+ using ThisType = EntropyBasedFluxEntropyCoordsFunction;
+
+public:
+ using GridViewType = GridViewImp;
+ using MomentBasis = MomentBasisImp;
+ using IndexSetType = typename GridViewType::IndexSet;
+ static const size_t dimFlux = MomentBasis::dimFlux;
+ static const size_t basis_dimRange = MomentBasis::dimRange;
+ using typename BaseType::DomainType;
+ using typename BaseType::E;
+ using typename BaseType::LocalFunctionType;
+ using typename BaseType::RangeFieldType;
+ using typename BaseType::StateType;
+ using ImplementationType = EntropyBasedFluxImplementation<MomentBasis>;
+ using AlphaReturnType = typename ImplementationType::AlphaReturnType;
+ using VectorType = typename ImplementationType::VectorType;
+ using I = XT::Grid::extract_intersection_t<GridViewType>;
+
+ explicit EntropyBasedFluxEntropyCoordsFunction(
+ const MomentBasis& basis_functions,
+ const RangeFieldType tau = 1e-9,
+ const RangeFieldType epsilon_gamma = 0.01,
+ const RangeFieldType chi = 0.5,
+ const RangeFieldType xi = 1e-3,
+ const std::vector<RangeFieldType> r_sequence = {0, 1e-8, 1e-6, 1e-4, 1e-3, 1e-2, 5e-2, 0.1, 0.5, 1},
+ const size_t k_0 = 500,
+ const size_t k_max = 1000,
+ const RangeFieldType epsilon = std::pow(2, -52))
+ : implementation_(std::make_shared<ImplementationType>(
+ basis_functions, tau, epsilon_gamma, chi, xi, r_sequence, k_0, k_max, epsilon))
+ {}
+
+ explicit EntropyBasedFluxEntropyCoordsFunction(EntropyBasedFluxFunction<GridViewType, MomentBasis>& other)
+ : implementation_(other.implementation_)
+ {}
+
+
+ static const constexpr bool available = true;
+
+ class Localfunction : public LocalFunctionType
+ {
+ using BaseType = LocalFunctionType;
+
+ public:
+ using typename BaseType::DynamicJacobianRangeType;
+ using typename BaseType::E;
+ using typename BaseType::RangeReturnType;
+
+ Localfunction(const ImplementationType& implementation)
+ : implementation_(implementation)
+ {}
+
+ virtual int order(const XT::Common::Parameter&) const override final
+ {
+ return 1.;
+ }
+
+ virtual RangeReturnType evaluate(const DomainType& /*point_in_reference_element*/,
+ const StateType& alpha,
+ const XT::Common::Parameter& /*param*/ = {}) const override final
+ {
+ return implementation_.evaluate_with_alpha(alpha);
+ }
+
+ virtual void jacobian(const DomainType& /*point_in_reference_element*/,
+ const StateType& alpha,
+ DynamicJacobianRangeType& result,
+ const XT::Common::Parameter& /*param*/ = {}) const override final
+ {
+ implementation_.jacobian_with_alpha(alpha, result);
+ } // ... jacobian(...)
+
+ private:
+ const ImplementationType& implementation_;
+ }; // class Localfunction
+
+ virtual bool x_dependent() const override final
+ {
+ return false;
+ }
+
+ virtual std::unique_ptr<LocalFunctionType> local_function() const override final
+ {
+ return std::make_unique<Localfunction>(*implementation_);
+ }
+
+ virtual std::unique_ptr<Localfunction> derived_local_function() const
+ {
+ return std::make_unique<Localfunction>(*implementation_);
+ }
+
+ StateType evaluate_kinetic_flux_precomputed(const StateType& flux_1,
+ const StateType& flux_2,
+ const DomainType& n_ij,
+ const size_t dd) const
+ {
+ return (flux_1 + flux_2) * n_ij[dd];
+ } // StateType evaluate_kinetic_flux(...)
+
+ StateType evaluate_kinetic_flux(const E& /*inside_entity*/,
+ const E& /*outside_entity*/,
+ const StateType& alpha_i,
+ const StateType& alpha_j,
+ const DomainType& n_ij,
+ const size_t dd) const
+ {
+ return implementation_->evaluate_kinetic_flux_with_alphas(alpha_i, alpha_j, n_ij, dd);
+ } // StateType evaluate_kinetic_flux(...)
+
+ const MomentBasis& basis_functions() const
+ {
+ return implementation_->basis_functions();
+ }
+
+ StateType get_u(const StateType& alpha) const
+ {
+ return implementation_->get_u(alpha);
+ }
+
+ StateType get_alpha(const StateType& u) const
+ {
+ const auto alpha = implementation_->get_alpha(u)->first;
+ StateType ret;
+ std::copy(alpha.begin(), alpha.end(), ret.begin());
+ return ret;
+ }
+
+ template <class FluxesMapType>
+ void calculate_minmod_density_reconstruction(const FieldVector<StateType, 3>& alphas,
+ FluxesMapType& precomputed_fluxes,
+ const size_t dd) const
+ {
+ implementation_->calculate_minmod_density_reconstruction(alphas, precomputed_fluxes, dd);
+ }
+
+ StateType calculate_boundary_flux(const StateType& alpha, const I& intersection)
+ {
+ return implementation_->calculate_boundary_flux(alpha, intersection);
+ }
+
+ void apply_inverse_hessian(const StateType& alpha, const StateType& u, StateType& Hinv_u) const
+ {
+ implementation_->apply_inverse_hessian(alpha, u, Hinv_u);
+ }
+
+private:
+ std::shared_ptr<ImplementationType> implementation_;
+};
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_LOCAL_FLUXES_ENTROPYBASED_KINETICCOORDS_HH
diff --git a/dune/gdt/momentmodels/hessianinverter.hh b/dune/gdt/momentmodels/hessianinverter.hh
new file mode 100644
index 0000000000000000000000000000000000000000..12f49ea7d406ba994a51b30dded5a4921ec347a6
--- /dev/null
+++ b/dune/gdt/momentmodels/hessianinverter.hh
@@ -0,0 +1,199 @@
+// 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:
+// Tobias Leibner (2018)
+
+#ifndef DUNE_GDT_MOMENTMODELS_HESSIANINVERTER_HH
+#define DUNE_GDT_MOMENTMODELS_HESSIANINVERTER_HH
+
+#include <string>
+
+#include <dune/xt/grid/functors/interfaces.hh>
+#include <dune/xt/common/parameter.hh>
+
+#include <dune/gdt/discretefunction/default.hh>
+#include <dune/gdt/momentmodels/entropyflux_kineticcoords.hh>
+#include <dune/gdt/operators/interfaces.hh>
+#include <dune/gdt/type_traits.hh>
+
+namespace Dune {
+namespace GDT {
+
+
+template <class SpaceType, class VectorType, class MomentBasis>
+class LocalEntropicHessianInverter : public XT::Grid::ElementFunctor<typename SpaceType::GridViewType>
+{
+ using GridViewType = typename SpaceType::GridViewType;
+ using BaseType = XT::Grid::ElementFunctor<GridViewType>;
+ using EntityType = typename GridViewType::template Codim<0>::Entity;
+ using IndexSetType = typename GridViewType::IndexSet;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<GridViewType, MomentBasis>;
+ using RangeFieldType = typename EntropyFluxType::RangeFieldType;
+ using LocalVectorType = typename EntropyFluxType::VectorType;
+ static const size_t dimFlux = EntropyFluxType::dimFlux;
+ static const size_t dimRange = EntropyFluxType::basis_dimRange;
+ using DiscreteFunctionType = DiscreteFunction<VectorType, GridViewType, dimRange, 1, RangeFieldType>;
+ using ConstDiscreteFunctionType = ConstDiscreteFunction<VectorType, GridViewType, dimRange, 1, RangeFieldType>;
+
+public:
+ explicit LocalEntropicHessianInverter(const SpaceType& space,
+ const VectorType& alpha_dofs,
+ const VectorType& u_update_dofs,
+ VectorType& alpha_range_dofs,
+ const EntropyFluxType& analytical_flux,
+ const RangeFieldType min_acceptable_density,
+ const XT::Common::Parameter& param,
+ const std::string filename = "")
+ : space_(space)
+ , alpha_(space_, alpha_dofs, "alpha")
+ , u_update_(space_, u_update_dofs, "u_update")
+ , local_alpha_(alpha_.local_discrete_function())
+ , local_u_update_(u_update_.local_discrete_function())
+ , range_(space_, alpha_range_dofs, "range")
+ , local_range_(range_.local_discrete_function())
+ , analytical_flux_(analytical_flux)
+ , local_flux_(analytical_flux_.derived_local_function())
+ , min_acceptable_density_(min_acceptable_density)
+ , param_(param)
+ , filename_(filename.empty() ? "" : (filename + "_regularization.txt"))
+ , index_set_(space_.grid_view().indexSet())
+ {}
+
+ explicit LocalEntropicHessianInverter(LocalEntropicHessianInverter& other)
+ : BaseType(other)
+ , space_(other.space_)
+ , alpha_(space_, other.alpha_.dofs().vector(), "source")
+ , u_update_(space_, other.u_update_.dofs().vector(), "source")
+ , local_alpha_(alpha_.local_discrete_function())
+ , local_u_update_(u_update_.local_discrete_function())
+ , range_(space_, other.range_.dofs().vector(), "range")
+ , local_range_(range_.local_discrete_function())
+ , analytical_flux_(other.analytical_flux_)
+ , local_flux_(analytical_flux_.derived_local_function())
+ , min_acceptable_density_(other.min_acceptable_density_)
+ , param_(other.param_)
+ , filename_(other.filename_)
+ , index_set_(space_.grid_view().indexSet())
+ {}
+
+ virtual XT::Grid::ElementFunctor<GridViewType>* copy() override final
+ {
+ return new LocalEntropicHessianInverter(*this);
+ }
+
+ void apply_local(const EntityType& entity) override final
+ {
+ local_alpha_->bind(entity);
+ local_u_update_->bind(entity);
+ local_range_->bind(entity);
+ XT::Common::FieldVector<RangeFieldType, dimRange> u, Hinv_u, alpha;
+ for (size_t ii = 0; ii < dimRange; ++ii) {
+ u[ii] = local_u_update_->dofs().get_entry(ii);
+ alpha[ii] = local_alpha_->dofs().get_entry(ii);
+ }
+ analytical_flux_.apply_inverse_hessian(alpha, u, Hinv_u);
+ for (auto&& entry : Hinv_u)
+ if (std::isnan(entry) || std::isinf(entry)) {
+ // std::cout << "x: " << entity.geometry().center() << "u: " << u << ", alpha: " << alpha << ", Hinv_u: "
+ // << Hinv_u << std::endl;
+ DUNE_THROW(Dune::MathError, "Hessian");
+ }
+
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ local_range_->dofs().set_entry(ii, Hinv_u[ii]);
+ } // void apply_local(...)
+
+private:
+ const SpaceType& space_;
+ const ConstDiscreteFunctionType alpha_;
+ const ConstDiscreteFunctionType u_update_;
+ std::unique_ptr<typename ConstDiscreteFunctionType::ConstLocalDiscreteFunctionType> local_alpha_;
+ std::unique_ptr<typename ConstDiscreteFunctionType::ConstLocalDiscreteFunctionType> local_u_update_;
+ DiscreteFunctionType range_;
+ std::unique_ptr<typename DiscreteFunctionType::LocalDiscreteFunctionType> local_range_;
+ const EntropyFluxType& analytical_flux_;
+ std::unique_ptr<typename EntropyFluxType::Localfunction> local_flux_;
+ const RangeFieldType min_acceptable_density_;
+ const XT::Common::Parameter& param_;
+ const std::string filename_;
+ const typename SpaceType::GridViewType::IndexSet& index_set_;
+}; // class LocalEntropicHessianInverter<...>
+
+template <class MomentBasisImp,
+ class SpaceImp,
+ class MatrixType = typename XT::LA::Container<typename MomentBasisImp::RangeFieldType>::MatrixType>
+class EntropicHessianInverter
+ : public OperatorInterface<MatrixType, typename SpaceImp::GridViewType, MomentBasisImp::dimRange, 1>
+{
+ using BaseType = OperatorInterface<MatrixType, typename SpaceImp::GridViewType, MomentBasisImp::dimRange, 1>;
+
+public:
+ using typename BaseType::VectorType;
+ using MomentBasis = MomentBasisImp;
+ using SpaceType = SpaceImp;
+ using SourceSpaceType = SpaceImp;
+ using RangeSpaceType = SpaceImp;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<typename SpaceType::GridViewType, MomentBasis>;
+ using RangeFieldType = typename MomentBasis::RangeFieldType;
+ using LocalVectorType = typename EntropyFluxType::VectorType;
+
+ EntropicHessianInverter(const EntropyFluxType& analytical_flux,
+ const SpaceType& space,
+ const RangeFieldType min_acceptable_density,
+ const std::string filename = "")
+ : analytical_flux_(analytical_flux)
+ , space_(space)
+ , min_acceptable_density_(min_acceptable_density)
+ , filename_(filename)
+ {}
+
+ virtual bool linear() const override final
+ {
+ return false;
+ }
+
+ virtual const SourceSpaceType& source_space() const override final
+ {
+ return space_;
+ }
+
+ virtual const RangeSpaceType& range_space() const override final
+ {
+ return space_;
+ }
+
+ void apply(const VectorType& /*source*/,
+ VectorType& /*range*/,
+ const XT::Common::Parameter& /*param*/) const override final
+ {
+ DUNE_THROW(Dune::NotImplemented, "Use apply_inverse_hessian!");
+ } // void apply(...)
+
+ void apply_inverse_hessian(const VectorType& alpha,
+ const VectorType& u_update,
+ VectorType& alpha_update,
+ const XT::Common::Parameter& param) const
+ {
+ LocalEntropicHessianInverter<SpaceType, VectorType, MomentBasis> local_hessian_inverter(
+ space_, alpha, u_update, alpha_update, analytical_flux_, min_acceptable_density_, param, filename_);
+ auto walker = XT::Grid::Walker<typename SpaceType::GridViewType>(space_.grid_view());
+ walker.append(local_hessian_inverter);
+ walker.walk(true);
+ } // void apply(...)
+
+private:
+ const EntropyFluxType& analytical_flux_;
+ const SpaceType& space_;
+ const RangeFieldType min_acceptable_density_;
+ const std::string filename_;
+}; // class EntropicHessianInverter<...>
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_MOMENTMODELS_HESSIANINVERTER_HH
diff --git a/dune/gdt/operators/advection-fv-entropybased.hh b/dune/gdt/operators/advection-fv-entropybased.hh
index 98392af3b1d0dcb8b2f5e076b86bcde65901b06a..e9a5b6fe1e09df3d279675dab5175645a3f52be0 100644
--- a/dune/gdt/operators/advection-fv-entropybased.hh
+++ b/dune/gdt/operators/advection-fv-entropybased.hh
@@ -17,6 +17,111 @@ namespace Dune {
namespace GDT {
+template <class OperatorImp, class InverseHessianOperatorImp>
+class EntropicCoordinatesOperator
+ : public OperatorInterface<typename OperatorImp::MatrixType, typename OperatorImp::SGV, OperatorImp::s_r>
+{
+ using BaseType = OperatorInterface<typename OperatorImp::MatrixType, typename OperatorImp::SGV, OperatorImp::s_r>;
+
+public:
+ using typename BaseType::RangeSpaceType;
+ using typename BaseType::SourceSpaceType;
+ using typename BaseType::VectorType;
+
+ using OperatorType = OperatorImp;
+ using InverseHessianOperatorType = InverseHessianOperatorImp;
+
+ EntropicCoordinatesOperator(const OperatorType& operator_in,
+ const InverseHessianOperatorType& inverse_hessian_operator)
+ : operator_(operator_in)
+ , inverse_hessian_operator_(inverse_hessian_operator)
+ {}
+
+ virtual bool linear() const override final
+ {
+ return false;
+ }
+
+ virtual const SourceSpaceType& source_space() const override final
+ {
+ return operator_.source_space();
+ }
+
+ virtual const RangeSpaceType& range_space() const override final
+ {
+ return operator_.range_space();
+ }
+
+ void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param) const override final
+ {
+ VectorType u_update = range;
+ std::fill(u_update.begin(), u_update.end(), 0.);
+ operator_.apply(source, u_update, param);
+ inverse_hessian_operator_.apply_inverse_hessian(source, u_update, range, param);
+ }
+
+ const OperatorType& operator_;
+ const InverseHessianOperatorType& inverse_hessian_operator_;
+}; // class EntropicCoordinatesOperator<...>
+
+template <class AdvectionOperatorImp, class RhsOperatorImp, class InverseHessianOperatorImp>
+class EntropicCoordinatesCombinedOperator
+ : public OperatorInterface<typename AdvectionOperatorImp::MatrixType,
+ typename AdvectionOperatorImp::SGV,
+ AdvectionOperatorImp::s_r>
+{
+ using BaseType = OperatorInterface<typename AdvectionOperatorImp::MatrixType,
+ typename AdvectionOperatorImp::SGV,
+ AdvectionOperatorImp::s_r>;
+
+public:
+ using typename BaseType::RangeSpaceType;
+ using typename BaseType::SourceSpaceType;
+ using typename BaseType::VectorType;
+
+ using AdvectionOperatorType = AdvectionOperatorImp;
+ using RhsOperatorType = RhsOperatorImp;
+ using InverseHessianOperatorType = InverseHessianOperatorImp;
+
+ EntropicCoordinatesCombinedOperator(const AdvectionOperatorType& advection_op,
+ const RhsOperatorType& rhs_op,
+ const InverseHessianOperatorType& inverse_hessian_operator)
+ : advection_op_(advection_op)
+ , rhs_op_(rhs_op)
+ , inverse_hessian_operator_(inverse_hessian_operator)
+ {}
+
+ virtual bool linear() const override final
+ {
+ return false;
+ }
+
+ virtual const SourceSpaceType& source_space() const override final
+ {
+ return advection_op_.source_space();
+ }
+
+ virtual const RangeSpaceType& range_space() const override final
+ {
+ return advection_op_.range_space();
+ }
+
+ void apply(const VectorType& source, VectorType& range, const XT::Common::Parameter& param) const override final
+ {
+ VectorType u_update = range;
+ std::fill(u_update.begin(), u_update.end(), 0.);
+ advection_op_.apply(source, u_update, param);
+ u_update *= -1.;
+ rhs_op_.apply(source, u_update, param);
+ inverse_hessian_operator_.apply_inverse_hessian(source, u_update, range, param);
+ }
+
+ const AdvectionOperatorType& advection_op_;
+ const RhsOperatorType& rhs_op_;
+ const InverseHessianOperatorType& inverse_hessian_operator_;
+}; // class EntropicCoordinatesOperator<...>
+
+
template <class AdvectionOperatorImp, class EntropySolverImp>
class EntropyBasedMomentFvOperator
: public OperatorInterface<typename AdvectionOperatorImp::MatrixType,
diff --git a/dune/gdt/operators/reconstruction/internal.hh b/dune/gdt/operators/reconstruction/internal.hh
index 2df45463dd9a330fda0ee4f40a56fa50d65ad630..850095c7c27f7f4bf32e943b678212879629bfd0 100644
--- a/dune/gdt/operators/reconstruction/internal.hh
+++ b/dune/gdt/operators/reconstruction/internal.hh
@@ -117,6 +117,55 @@ constexpr size_t EigenvectorWrapperBase<AnalyticalFluxType, MatrixImp, VectorImp
template <class AnalyticalFluxType, class MatrixImp, class VectorImp>
constexpr size_t EigenvectorWrapperBase<AnalyticalFluxType, MatrixImp, VectorImp>::dimRange;
+// This class does not perform any computations, use this class if you want to reconstruct in ordinary coordinates
+// instead of characteristic coordinates
+template <class AnalyticalFluxType, class VectorImp>
+class DummyEigenVectorWrapper : public EigenvectorWrapperBase<AnalyticalFluxType, int, VectorImp>
+{
+ using BaseType = EigenvectorWrapperBase<AnalyticalFluxType, int, VectorImp>;
+
+public:
+ using typename BaseType::DomainType;
+ using typename BaseType::E;
+ using typename BaseType::MatrixType;
+ using typename BaseType::VectorType;
+
+ DummyEigenVectorWrapper(const AnalyticalFluxType& analytical_flux, const bool flux_is_affine)
+ : BaseType(analytical_flux, flux_is_affine)
+ {}
+
+ virtual void apply_eigenvectors(const size_t /*dd*/, const VectorType& u, VectorType& ret) const override final
+ {
+ ret = u;
+ }
+ virtual void
+ apply_inverse_eigenvectors(const size_t /*dd*/, const VectorType& u, VectorType& ret) const override final
+ {
+ ret = u;
+ }
+
+ virtual void compute_eigenvectors(const E& /*entity*/,
+ const DomainType& /*x_local*/,
+ const VectorType& /*u*/,
+ const XT::Common::Parameter& /*param*/) const override final
+ {}
+
+ virtual void compute_eigenvectors_impl(const E& /*entity*/,
+ const DomainType& /*x_local*/,
+ const VectorType& /*u*/,
+ const XT::Common::Parameter& /*param*/) const override final
+ {}
+
+ virtual const MatrixType& eigenvectors(const size_t /*dd*/) const override final
+ {
+ DUNE_THROW(Dune::NotImplemented, "This class does not provide eigenvectors!");
+ return zero_;
+ }
+
+private:
+ static const MatrixType zero_ = 0.;
+};
+
template <class AnalyticalFluxType,
class MatrixType = XT::LA::CommonDenseMatrix<typename AnalyticalFluxType::R>,
class VectorType = FieldVector<typename AnalyticalFluxType::RangeFieldType, AnalyticalFluxType::rC>>
diff --git a/dune/gdt/operators/reconstruction/linear.hh b/dune/gdt/operators/reconstruction/linear.hh
index 375214d16ad36f5fcb2236124d5350547ae3aea6..580d34b819df7fe673077bbbed8ce96a83c48b73 100644
--- a/dune/gdt/operators/reconstruction/linear.hh
+++ b/dune/gdt/operators/reconstruction/linear.hh
@@ -127,7 +127,7 @@ private:
const size_t dd = intersection.indexInInside() / 2;
const size_t index = (intersection.indexInInside() % 2) * 2;
if (intersection.boundary() && !intersection.neighbor()) // boundary intersections
- stencils_[dd][index] = boundary_values_.evaluate(entity.geometry().local(intersection.geometry().center()));
+ stencils_[dd][index] = boundary_values_.evaluate(intersection.geometry().center());
else if (intersection.neighbor()) // inner and periodic intersections
stencils_[dd][index] = source_values_[grid_view_.indexSet().index(intersection.outside())];
else if (!intersection.neighbor() && !intersection.boundary()) // processor boundary
diff --git a/dune/gdt/operators/reconstruction/linear_kinetic.hh b/dune/gdt/operators/reconstruction/linear_kinetic.hh
new file mode 100644
index 0000000000000000000000000000000000000000..30ec756c78498ad4d46daedb1e154ffbee83ffad
--- /dev/null
+++ b/dune/gdt/operators/reconstruction/linear_kinetic.hh
@@ -0,0 +1,190 @@
+// 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:
+// Tobias Leibner (2019)
+
+#ifndef DUNE_GDT_OPERATORS_FV_RECONSTRUCTION_LINEAR_KINETIC_HH
+#define DUNE_GDT_OPERATORS_FV_RECONSTRUCTION_LINEAR_KINETIC_HH
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/xt/common/fvector.hh>
+#include <dune/xt/common/parameter.hh>
+
+#include <dune/xt/grid/walker.hh>
+
+#include <dune/gdt/operators/interfaces.hh>
+#include <dune/gdt/spaces/l2/discontinuous-lagrange.hh>
+#include <dune/gdt/tools/discretevalued-grid-function.hh>
+
+#include "slopes.hh"
+#include "internal.hh"
+
+namespace Dune {
+namespace GDT {
+
+
+template <class GV, class AnalyticalFluxType, class BoundaryValueType, class LocalVectorType>
+class LocalPointwiseLinearKineticReconstructionOperator : public XT::Grid::ElementFunctor<GV>
+{
+ using ThisType = LocalPointwiseLinearKineticReconstructionOperator;
+ using BaseType = XT::Grid::ElementFunctor<GV>;
+ static constexpr size_t dimDomain = BoundaryValueType::d;
+ static constexpr size_t dimRange = BoundaryValueType::r;
+ using EntityType = typename GV::template Codim<0>::Entity;
+ using RangeFieldType = typename BoundaryValueType::RangeFieldType;
+ using StencilType = DynamicVector<LocalVectorType>;
+ using StencilsType = std::vector<StencilType>;
+ using DomainType = typename BoundaryValueType::DomainType;
+ using RangeType = typename BoundaryValueType::RangeReturnType;
+ static constexpr size_t stencil_size = 3;
+ using ReconstructedFunctionType = DiscreteValuedGridFunction<GV, dimRange, 1, RangeFieldType>;
+
+public:
+ explicit LocalPointwiseLinearKineticReconstructionOperator(ReconstructedFunctionType& reconstructed_function,
+ const GV& grid_view,
+ const AnalyticalFluxType& analytical_flux,
+ const std::vector<LocalVectorType>& source_values,
+ const BoundaryValueType& boundary_values,
+ const XT::Common::Parameter& param)
+ : reconstructed_function_(reconstructed_function)
+ , grid_view_(grid_view)
+ , analytical_flux_(analytical_flux)
+ , source_values_(source_values)
+ , boundary_values_(boundary_values)
+ , param_(param)
+ , stencils_(dimDomain, StencilType(stencil_size))
+ {}
+
+ LocalPointwiseLinearKineticReconstructionOperator(const ThisType& other)
+ : BaseType(other)
+ , reconstructed_function_(other.reconstructed_function_)
+ , grid_view_(other.grid_view_)
+ , analytical_flux_(other.analytical_flux_)
+ , source_values_(other.source_values_)
+ , boundary_values_(other.boundary_values_)
+ , param_(other.param_)
+ , stencils_(dimDomain, StencilType(stencil_size))
+ {}
+
+ virtual XT::Grid::ElementFunctor<GV>* copy() override final
+ {
+ return new ThisType(*this);
+ }
+
+ virtual void apply_local(const EntityType& entity) override final
+ {
+ // In a MPI parallel run, if entity is on boundary of overlap, we do not have to reconstruct
+ if (!fill_stencils(entity))
+ return;
+
+ auto& local_reconstructed_values = reconstructed_function_.local_values(entity);
+ for (size_t dd = 0; dd < dimDomain; ++dd)
+ analytical_flux_.calculate_minmod_density_reconstruction(stencils_[dd], local_reconstructed_values, dd);
+ } // void apply_local(...)
+
+private:
+ bool fill_stencils(const EntityType& entity)
+ {
+ const auto entity_index = grid_view_.indexSet().index(entity);
+ for (size_t dd = 0; dd < dimDomain; ++dd)
+ stencils_[dd][1] = source_values_[entity_index];
+ for (const auto& intersection : Dune::intersections(grid_view_, entity)) {
+ const size_t dd = intersection.indexInInside() / 2;
+ const size_t index = (intersection.indexInInside() % 2) * 2;
+ if (intersection.boundary() && !intersection.neighbor()) { // boundary intersections
+ stencils_[dd][index] = boundary_values_.evaluate(intersection.geometry().center());
+ } else if (intersection.neighbor()) { // inner and periodic intersections {
+ stencils_[dd][index] = source_values_[grid_view_.indexSet().index(intersection.outside())];
+ } else if (!intersection.neighbor() && !intersection.boundary()) { // processor boundary
+ return false;
+ } else {
+ DUNE_THROW(Dune::NotImplemented, "This should not happen!");
+ }
+ } // intersections
+ return true;
+ } // void fill_stencils(...)
+
+ ReconstructedFunctionType& reconstructed_function_;
+ const GV& grid_view_;
+ const AnalyticalFluxType& analytical_flux_;
+ const std::vector<LocalVectorType>& source_values_;
+ const BoundaryValueType& boundary_values_;
+ const XT::Common::Parameter& param_;
+ StencilsType stencils_;
+}; // class LocalPointwiseLinearKineticReconstructionOperator
+
+// Does not reconstruct a full first-order DG function, but only stores the reconstructed values at the intersection
+// centers. This avoids the interpolation in this operator and the evaluation of the reconstructed function in the
+// finite volume operator which are both quite expensive for large dimRange.
+template <class GV, class BoundaryValueImp, class AnalyticalFluxType, class VectorType, class LocalVectorType>
+class PointwiseLinearKineticReconstructionOperator
+{
+public:
+ using BoundaryValueType = BoundaryValueImp;
+ using E = XT::Grid::extract_entity_t<GV>;
+ using EigenVectorWrapperType = DummyEigenVectorWrapper<LocalVectorType>;
+ using R = typename BoundaryValueType::R;
+ static constexpr size_t dimDomain = BoundaryValueType::d;
+ static constexpr size_t dimRange = BoundaryValueType::r;
+ using SpaceType = SpaceInterface<GV, dimRange, 1, R>;
+ using ReconstructedFunctionType = DiscreteValuedGridFunction<GV, dimRange, 1, R>;
+ using ReconstructedValuesType = std::vector<typename ReconstructedFunctionType::LocalFunctionValuesType>;
+
+ PointwiseLinearKineticReconstructionOperator(const BoundaryValueType& boundary_values,
+ const SpaceType& space,
+ const AnalyticalFluxType& analytical_flux)
+ : boundary_values_(boundary_values)
+ , space_(space)
+ , analytical_flux_(analytical_flux)
+ {}
+
+ bool linear() const
+ {
+ return false;
+ }
+
+ const SpaceType& space() const
+ {
+ return space_;
+ }
+
+ void apply(const VectorType& source, ReconstructedFunctionType& range, const XT::Common::Parameter& param) const
+ {
+ // evaluate cell averages
+ const auto& grid_view = space_.grid_view();
+ const auto& index_set = grid_view.indexSet();
+ std::vector<LocalVectorType> source_values(index_set.size(0));
+ // Contains the evaluations of exp(alpha^T b(v_i)) at each quadrature point v_i
+ auto source_func = make_discrete_function(space_, source);
+ const auto local_source = source_func.local_function();
+ for (const auto& entity : Dune::elements(grid_view)) {
+ local_source->bind(entity);
+ const auto entity_index = index_set.index(entity);
+ source_values[entity_index] = local_source->evaluate(entity.geometry().local(entity.geometry().center()));
+ }
+
+ // do reconstruction
+ auto local_reconstruction_operator =
+ LocalPointwiseLinearKineticReconstructionOperator<GV, AnalyticalFluxType, BoundaryValueType, LocalVectorType>(
+ range, grid_view, analytical_flux_, source_values, boundary_values_, param);
+ auto walker = XT::Grid::Walker<GV>(grid_view);
+ walker.append(local_reconstruction_operator);
+ walker.walk(true);
+ } // void apply(...)
+
+private:
+ const BoundaryValueType& boundary_values_;
+ const SpaceType& space_;
+ const AnalyticalFluxType& analytical_flux_;
+}; // class PointwiseLinearKineticReconstructionOperator<...>
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_OPERATORS_FV_RECONSTRUCTION_LINEAR_KINETIC_HH
diff --git a/dune/gdt/operators/reconstruction/slopes.hh b/dune/gdt/operators/reconstruction/slopes.hh
index 6a35099e8e4ed39dec396107e418dc940cfe38fd..a840a2aa7afe8fb16fa298255037b8cdeba68728 100644
--- a/dune/gdt/operators/reconstruction/slopes.hh
+++ b/dune/gdt/operators/reconstruction/slopes.hh
@@ -213,6 +213,16 @@ public:
return minmod(slope_left_char, slope_right_char);
}
+ virtual VectorType get(const E& /*entity*/,
+ const StencilType& stencil,
+ const EigenVectorWrapperType& /*eigenvectors*/,
+ const size_t /*dd*/) const override final
+ {
+ const VectorType slope_left = stencil[1] - stencil[0];
+ const VectorType slope_right = stencil[2] - stencil[1];
+ return minmod(slope_left, slope_right);
+ }
+
static VectorType minmod(const VectorType& first_slope, const VectorType& second_slope)
{
VectorType ret(0.);
diff --git a/dune/gdt/test/entropic-coords-mn-discretization.hh b/dune/gdt/test/entropic-coords-mn-discretization.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b6b24fa666dffd0dd524f3a73f195339bf0efafc
--- /dev/null
+++ b/dune/gdt/test/entropic-coords-mn-discretization.hh
@@ -0,0 +1,384 @@
+// This file is part of the dune-gdt project:
+// https://github.com/dune-community/dune-gdt
+// Copyright 2010-2016 dune-gdt developers and contributors. All rights reserved.
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// Authors:
+// Tobias Leibner (2016)
+
+#ifndef DUNE_GDT_TEST_HYPERBOLIC_ENTROPIC_COORDS_MN_DISCRETIZATION_HH
+#define DUNE_GDT_TEST_HYPERBOLIC_ENTROPIC_COORDS_MN_DISCRETIZATION_HH
+
+#include <chrono>
+
+#include <dune/xt/common/parallel/threadmanager.hh>
+#include <dune/xt/common/string.hh>
+#include <dune/xt/common/test/gtest/gtest.h>
+
+#include <dune/xt/grid/information.hh>
+#include <dune/xt/grid/gridprovider.hh>
+
+#include <dune/xt/la/container.hh>
+
+#include <dune/gdt/discretefunction/default.hh>
+#include <dune/gdt/operators/advection-fv-entropybased.hh>
+#include <dune/gdt/operators/advection-fv.hh>
+#include <dune/gdt/operators/generic.hh>
+#include <dune/gdt/operators/reconstruction/linear_nocharacteristiccoords.hh>
+#include <dune/gdt/operators/reconstruction/linear_kinetic.hh>
+#include <dune/gdt/interpolations/default.hh>
+#include <dune/gdt/momentmodels/hessianinverter.hh>
+#include <dune/gdt/momentmodels/entropyflux_kineticcoords.hh>
+#include <dune/gdt/momentmodels/entropyflux.hh>
+#include <dune/gdt/local/numerical-fluxes/kinetic.hh>
+#include <dune/gdt/local/operators/advection-fv.hh>
+#include <dune/gdt/spaces/l2/finite-volume.hh>
+#include <dune/gdt/tools/timestepper/fractional-step.hh>
+#include <dune/gdt/tools/timestepper/explicit-rungekutta.hh>
+#include <dune/gdt/tools/timestepper/adaptive-rungekutta.hh>
+#include <dune/gdt/tools/timestepper/matrix-exponential-kinetic-isotropic.hh>
+
+#include <dune/gdt/test/momentmodels/kineticequation.hh>
+
+#include "pn-discretization.hh"
+
+template <class TestCaseType>
+struct HyperbolicEntropicCoordsMnDiscretization
+{
+ // returns: (l1norm, l2norm, linfnorm, MPI rank)
+ static std::pair<Dune::FieldVector<double, 3>, int> run(size_t num_save_steps = 1,
+ size_t num_output_steps = 0,
+ size_t quad_order = size_t(-1),
+ size_t quad_refinements = size_t(-1),
+ std::string grid_size = "",
+ size_t overlap_size = 2,
+ double t_end = 0.,
+ std::string filename = "",
+ bool /*disable_thread_cache*/ = false)
+ {
+ using namespace Dune;
+ using namespace Dune::GDT;
+
+ //******************* get typedefs and constants from ProblemType **********************//
+ using MomentBasis = typename TestCaseType::MomentBasis;
+ using DiscreteFunctionType = typename TestCaseType::DiscreteFunctionType;
+ using GridType = typename TestCaseType::GridType;
+ using SpaceType = typename TestCaseType::SpaceType;
+ using AdvectionSourceSpaceType = typename TestCaseType::AdvectionSourceSpaceType;
+ using GV = typename TestCaseType::GridViewType;
+ // using E = XT::Grid::extract_entity_t<GV>;
+ using I = XT::Grid::extract_intersection_t<GV>;
+ using ProblemType = typename TestCaseType::ProblemType;
+ using RangeFieldType = typename MomentBasis::RangeFieldType;
+ using BoundaryValueType = typename ProblemType::BoundaryValueType;
+ static constexpr size_t dimDomain = MomentBasis::dimDomain;
+ static constexpr size_t dimRange = MomentBasis::dimRange;
+ using MatrixType = typename XT::LA::Container<RangeFieldType>::MatrixType;
+ using VectorType = typename XT::LA::Container<RangeFieldType>::VectorType;
+ using GenericFunctionType = XT::Functions::GenericFunction<dimDomain, dimRange, 1, RangeFieldType>;
+ using DomainType = FieldVector<RangeFieldType, dimDomain>;
+ using RangeType = FieldVector<RangeFieldType, dimRange>;
+ // static const RangeFieldType scale_factor = 1e2;
+
+ //******************* create grid and FV space ***************************************
+ auto grid_config = ProblemType::default_grid_cfg();
+ if (!grid_size.empty())
+ grid_config["num_elements"] = grid_size;
+ grid_config["overlap_size"] = XT::Common::to_string(overlap_size);
+ const auto grid_ptr =
+ Dune::XT::Grid::CubeGridProviderFactory<GridType>::create(grid_config, MPIHelper::getCommunicator()).grid_ptr();
+ assert(grid_ptr->comm().size() == 1 || grid_ptr->overlapSize(0) > 0);
+ const GV grid_view(grid_ptr->leafGridView());
+ const SpaceType fv_space(grid_view);
+ const AdvectionSourceSpaceType advection_source_space(grid_view, 1);
+
+ //******************* create EquationType object ***************************************
+ std::shared_ptr<const MomentBasis> basis_functions = std::make_shared<const MomentBasis>(
+ quad_order == size_t(-1) ? MomentBasis::default_quad_order() : quad_order,
+ quad_refinements == size_t(-1) ? MomentBasis::default_quad_refinements() : quad_refinements);
+ const std::unique_ptr<ProblemType> problem_ptr =
+ XT::Common::make_unique<ProblemType>(*basis_functions, grid_view, grid_config);
+ const auto& problem = *problem_ptr;
+ const auto initial_values_u = problem.initial_values();
+ const auto boundary_values_u = problem.boundary_values();
+
+ using AnalyticalFluxType = typename ProblemType::FluxType;
+ using EntropyFluxType = EntropyBasedFluxEntropyCoordsFunction<GV, MomentBasis>;
+ using OldEntropyFluxType = EntropyBasedFluxFunction<GV, MomentBasis>;
+ 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
+ std::map<DomainType, RangeType, XT::Common::FieldVectorFloatLess> alpha_boundary_vals;
+ for (const auto& element : Dune::elements(grid_view))
+ for (const auto& intersection : Dune::intersections(grid_view, element))
+ if (intersection.boundary()) {
+ const auto x = intersection.geometry().center();
+ const auto u = boundary_values_u->evaluate(x);
+ alpha_boundary_vals.insert(std::make_pair(x, analytical_flux->get_alpha(u)));
+ }
+ GenericFunctionType boundary_values_alpha(
+ 1, [&](const DomainType& x, const XT::Common::Parameter&) { return alpha_boundary_vals[x]; });
+
+ // calculate boundary kinetic fluxes
+ std::map<DomainType, RangeType, XT::Common::FieldVectorFloatLess> boundary_fluxes;
+ for (const auto& element : Dune::elements(grid_view))
+ for (const auto& intersection : Dune::intersections(grid_view, element))
+ if (intersection.boundary()) {
+ const auto x = intersection.geometry().center();
+ const auto boundary_flux = analytical_flux->calculate_boundary_flux(alpha_boundary_vals[x], intersection);
+ boundary_fluxes.insert(std::make_pair(x, boundary_flux));
+ }
+ GenericFunctionType boundary_kinetic_fluxes(
+ 1, [&](const DomainType& x, const XT::Common::Parameter&) { return boundary_fluxes[x]; });
+
+
+ // ***************** project initial values to discrete function *********************
+ // create a discrete function for the solution
+ DiscreteFunctionType u(fv_space, "solution");
+ DiscreteFunctionType alpha(fv_space, "solution");
+ // project initial values
+ default_interpolation(*initial_values_u, u, grid_view);
+
+ // convert initial values to alpha
+ 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;
+ for (auto&& element : Dune::elements(grid_view)) {
+ u_local_func->bind(element);
+ alpha_local_func->bind(element);
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ u_local[ii] = u_local_func->dofs().get_entry(ii);
+ const auto alpha_local = analytical_flux->get_alpha(u_local);
+ for (auto&& entry : alpha_local)
+ if (entry > 1) {
+ std::cout << XT::Common::to_string(element.geometry().center()) << ", " << XT::Common::to_string(u_local)
+ << ", " << XT::Common::to_string(alpha_local) << std::endl;
+ break;
+ }
+ for (size_t ii = 0; ii < dimRange; ++ii)
+ alpha_local_func->dofs().set_entry(ii, alpha_local[ii]);
+ }
+
+ // ******************** choose flux and rhs operator and timestepper ******************************************
+
+ using AdvectionOperatorType = AdvectionFvOperator<MatrixType, GV, dimRange>;
+ using HessianInverterType = EntropicHessianInverter<MomentBasis, SpaceType>;
+#if 0
+ using ReconstructionOperatorType = PointwiseLinearReconstructionNoCharOperator<
+ GV,
+ BoundaryValueType,
+ EntropyFluxType,
+ VectorType,
+ RangeType>;
+#else
+ using ReconstructionOperatorType =
+ PointwiseLinearKineticReconstructionOperator<GV, BoundaryValueType, EntropyFluxType, VectorType, RangeType>;
+#endif
+ using ReconstructionAdvectionOperatorType =
+ AdvectionWithPointwiseReconstructionOperator<AdvectionOperatorType, ReconstructionOperatorType>;
+ using NonEntropicAdvectionOperatorType =
+ std::conditional_t<TestCaseType::reconstruction, ReconstructionAdvectionOperatorType, AdvectionOperatorType>;
+ // using FvOperatorType = EntropicCoordinatesOperator<
+ // NonEntropicAdvectionOperatorType,
+ // HessianInverterType>;
+ using NonEntropicRhsOperatorType = GenericOperator<MatrixType, GV, dimRange>;
+ // using RhsOperatorType = EntropicCoordinatesOperator<NonEntropicRhsOperatorType, HessianInverterType>;
+ using CombinedOperatorType = EntropicCoordinatesCombinedOperator<NonEntropicAdvectionOperatorType,
+ NonEntropicRhsOperatorType,
+ HessianInverterType>;
+
+ // using OperatorTimeStepperType =
+ // ExplicitRungeKuttaTimeStepper<FvOperatorType,
+ // DiscreteFunctionType,
+ // TimeStepperMethods::explicit_euler>;
+ // using RhsTimeStepperType =
+ // ExplicitRungeKuttaTimeStepper<RhsOperatorType,
+ // DiscreteFunctionType,
+ // TimeStepperMethods::explicit_euler>;
+
+ // using OperatorTimeStepperType =
+ // AdaptiveRungeKuttaTimeStepper<FvOperatorType,
+ // DiscreteFunctionType>;
+ // using RhsTimeStepperType =
+ // AdaptiveRungeKuttaTimeStepper<RhsOperatorType,
+ // DiscreteFunctionType>;
+
+ // using TimeStepperType = FractionalTimeStepper<OperatorTimeStepperType, RhsTimeStepperType>;
+
+ using TimeStepperType =
+ AdaptiveRungeKuttaTimeStepper<CombinedOperatorType, DiscreteFunctionType, TimeStepperMethods::dormand_prince>;
+
+ // *************** Calculate dx and initial dt **************************************
+ Dune::XT::Grid::Dimensions<GV> dimensions(grid_view);
+ RangeFieldType dx = dimensions.entity_width.max();
+ if (dimDomain == 2)
+ dx /= std::sqrt(2);
+ if (dimDomain == 3)
+ dx /= std::sqrt(3);
+ RangeFieldType dt = CFL * dx;
+
+ // *********************** create operators and timesteppers ************************************
+ NumericalKineticFlux<GV, MomentBasis, EntropyFluxType> numerical_flux(*analytical_flux, *basis_functions);
+ AdvectionOperatorType advection_operator(grid_view, numerical_flux, advection_source_space, fv_space);
+
+ // boundary treatment
+ using BoundaryOperator =
+ LocalAdvectionFvBoundaryTreatmentByCustomExtrapolationOperator<I, VectorType, GV, dimRange>;
+ using LambdaType = typename BoundaryOperator::LambdaType;
+ using StateType = typename BoundaryOperator::StateType;
+
+ LambdaType boundary_lambda =
+ [&](const I& intersection,
+ const FieldVector<RangeFieldType, dimDomain - 1>& xx_in_reference_intersection_coordinates,
+ const AnalyticalFluxType& /*flux*/,
+ const StateType& /*u*/,
+ const XT::Common::Parameter& /*param*/) {
+ // return
+ // boundary_values_alpha.evaluate(intersection.geometry().global(xx_in_reference_intersection_coordinates));
+ return boundary_kinetic_fluxes.evaluate(
+ intersection.geometry().global(xx_in_reference_intersection_coordinates));
+ };
+ XT::Grid::ApplyOn::NonPeriodicBoundaryIntersections<GV> filter;
+ advection_operator.append(boundary_lambda, {}, filter);
+
+ // constexpr double epsilon = 1e-11;
+ // MinmodSlope<E, DummyEigenVectorWrapper<RangeType>> slope;
+ // ReconstructionOperatorType reconstruction_operator(boundary_values_alpha, fv_space, *analytical_flux, slope);
+ ReconstructionOperatorType reconstruction_operator(boundary_values_alpha, fv_space, *analytical_flux);
+ ReconstructionAdvectionOperatorType reconstruction_advection_operator(advection_operator, reconstruction_operator);
+
+ if (XT::Common::is_zero(t_end))
+ t_end = problem.t_end();
+
+ if (!filename.empty())
+ filename += "_minmod_";
+ filename += ProblemType::static_id();
+ filename += "_grid_" + grid_config["num_elements"];
+ filename += "_tend_" + XT::Common::to_string(t_end);
+ filename += "_quad_" + XT::Common::to_string(quad_order);
+ filename += TestCaseType::reconstruction ? "_ord2" : "_ord1";
+ filename += "_" + basis_functions->mn_name();
+
+ HessianInverterType hessian_inverter(
+ *analytical_flux, fv_space, problem.psi_vac() * basis_functions->unit_ball_volume() / 10, filename);
+ auto& non_entropic_advection_operator =
+ FvOperatorChooser<TestCaseType::reconstruction>::choose(advection_operator, reconstruction_advection_operator);
+ // FvOperatorType fv_operator(non_entropic_advection_operator, hessian_inverter);
+
+
+ const auto u_iso = basis_functions->u_iso();
+ const auto basis_integrated = basis_functions->integrated();
+ const auto sigma_a = problem.sigma_a();
+ const auto sigma_s = problem.sigma_s();
+ const auto Q = problem.Q();
+ auto rhs_func = [&](const auto& /*source*/,
+ const auto& local_source,
+ auto& local_range,
+ const Dune::XT::Common::Parameter& /*param*/) {
+ const auto& element = local_range.element();
+ local_source->bind(element);
+ const auto center = element.geometry().center();
+ const auto alpha = local_source->evaluate(center);
+ const auto u = analytical_flux->get_u(alpha);
+ const auto sigma_a_value = sigma_a->evaluate(center)[0];
+ const auto sigma_s_value = sigma_s->evaluate(center)[0];
+ const auto sigma_t_value = sigma_a_value + sigma_s_value;
+ const auto Q_value = Q->evaluate(center)[0];
+ auto ret = u * (-sigma_t_value);
+ ret += u_iso * basis_functions->density(u) * sigma_s_value;
+ ret += basis_integrated * Q_value;
+ for (size_t ii = 0; ii < local_range.dofs().size(); ++ii)
+ local_range.dofs()[ii] += ret[ii];
+ };
+ NonEntropicRhsOperatorType non_entropic_rhs_operator(
+ fv_space, fv_space, std::vector<typename NonEntropicRhsOperatorType::GenericElementFunctionType>(1, rhs_func));
+ // RhsOperatorType rhs_operator(non_entropic_rhs_operator, hessian_inverter);
+ CombinedOperatorType combined_operator(
+ non_entropic_advection_operator, non_entropic_rhs_operator, hessian_inverter);
+
+ // ******************************** do the time steps ***********************************************************
+ // OperatorTimeStepperType timestepper_op(fv_operator, alpha, -1.0);
+ // RhsTimeStepperType timestepper_rhs(rhs_operator, alpha, 1.0);
+ // TimeStepperType timestepper(timestepper_op, timestepper_rhs);
+ TimeStepperType timestepper(combined_operator, alpha, 1.);
+
+ auto begin_time = std::chrono::steady_clock::now();
+ auto visualizer = std::make_unique<XT::Functions::GenericVisualizer<dimRange, 1, double>>(
+ 1, [&basis_functions, &analytical_flux](const int /*comp*/, const auto& val) {
+ return basis_functions->density(analytical_flux->get_u(val));
+ });
+ timestepper.solve(t_end,
+ dt,
+ num_save_steps,
+ num_output_steps,
+ false,
+ true,
+ false,
+ true,
+ false,
+ filename,
+ // *basis_functions->visualizer(),
+ *visualizer,
+ basis_functions->stringifier());
+ auto end_time = std::chrono::steady_clock::now();
+ std::chrono::duration<double> time_diff = end_time - begin_time;
+ if (grid_view.comm().rank() == 0)
+ std::cout << "Solving took: " << XT::Common::to_string(time_diff.count(), 15) << " s" << std::endl;
+
+ auto ret = std::make_pair(FieldVector<double, 3>(0.), int(0));
+ double& l1norm = ret.first[0];
+ double& l2norm = ret.first[1];
+ double& linfnorm = ret.first[2];
+ ret.second = grid_view.comm().rank();
+ const auto& current_sol = timestepper.current_solution();
+ const auto local_sol = current_sol.local_function();
+ for (const auto& entity : elements(grid_view, Dune::Partitions::interior)) {
+ local_sol->bind(entity);
+ const auto val = local_sol->evaluate(entity.geometry().local(entity.geometry().center()));
+ RangeFieldType psi = basis_functions->density(val);
+ l1norm += std::abs(psi) * entity.geometry().volume();
+ l2norm += std::pow(psi, 2) * entity.geometry().volume();
+ linfnorm = std::max(std::abs(psi), linfnorm);
+ }
+ l1norm = grid_view.comm().sum(l1norm);
+ l2norm = grid_view.comm().sum(l2norm);
+ linfnorm = grid_view.comm().max(linfnorm);
+ l2norm = std::sqrt(l2norm);
+ return ret;
+ }
+};
+
+template <class TestCaseType>
+struct HyperbolicEntropicCoordsMnTest
+ : public HyperbolicEntropicCoordsMnDiscretization<TestCaseType>
+ , public ::testing::Test
+{
+ void run()
+ {
+ auto norms = HyperbolicEntropicCoordsMnDiscretization<TestCaseType>::run(
+ 10,
+ -1,
+ TestCaseType::RealizabilityLimiterChooserType::quad_order,
+ TestCaseType::RealizabilityLimiterChooserType::quad_refinements,
+ "",
+ 2,
+ TestCaseType::t_end,
+ "test_dormandprince",
+ Dune::GDT::is_full_moment_basis<typename TestCaseType::MomentBasis>::value)
+ .first;
+ const double l1norm = norms[0];
+ const double l2norm = norms[1];
+ const double linfnorm = norms[2];
+ using ResultsType = typename TestCaseType::ExpectedResultsType;
+ EXPECT_NEAR(ResultsType::l1norm, l1norm, ResultsType::l1norm * ResultsType::tol);
+ EXPECT_NEAR(ResultsType::l2norm, l2norm, ResultsType::l2norm * ResultsType::tol);
+ EXPECT_NEAR(ResultsType::linfnorm, linfnorm, ResultsType::linfnorm * ResultsType::tol);
+ }
+};
+
+#endif // DUNE_GDT_TEST_HYPERBOLIC_ENTROPIC_COORDS_MN_DISCRETIZATION_HH
diff --git a/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc b/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c805ad1656256b13e21cc6d87a8ccc5a2beb7b7a
--- /dev/null
+++ b/dune/gdt/test/hyperbolic__momentmodels__entropic_coords_mn.cc
@@ -0,0 +1,63 @@
+// This file is part of the dune-gdt project:
+// https://github.com/dune-community/dune-gdt
+// Copyright 2010-2016 dune-gdt developers and contributors. All rights reserved.
+// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+// Authors:
+// Tobias Leibner (2016)
+
+// This one has to come first (includes the config.h)!
+#include <dune/xt/common/test/main.hxx>
+
+#define USE_LP_POSITIVITY_LIMITER 1
+
+#include <dune/gdt/test/momentmodels/kinetictransport/testcases.hh>
+#include <dune/gdt/test/entropic-coords-mn-discretization.hh>
+
+using Yasp1 = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>;
+using Yasp2 = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double, 2>>;
+using Yasp3 = Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double, 3>>;
+
+using YaspGridTestCasesAll = testing::Types<
+#if HAVE_CLP
+// Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 7>, false>,
+// Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 20>, true>
+// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 20>, true>
+// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::LegendreMomentBasis<double, double, 21>, true>
+#endif
+ // Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 8, 1, 1>, false>
+ Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 100, 1, 1>, true>
+// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 20, 1, 1>, false>,
+// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::HatFunctionMomentBasis<double, 1, double, 100, 1, 1>, true>
+// Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, false>
+// Dune::GDT::SourceBeamMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, true>,
+// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, false>
+// Dune::GDT::PlaneSourceMnTestCase<Yasp1, Dune::GDT::PartialMomentBasis<double, 1, double, 8, 1, 1>, true>
+#if !DXT_DISABLE_LARGE_TESTS
+// ,
+# if HAVE_CLP
+// Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2, 3>,
+// false>, Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2,
+// 3>, true> Dune::GDT::CheckerboardMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3>,
+// true> Dune::GDT::CheckerboardMnTestCase<Yasp3, Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2, 3>,
+// false>, Dune::GDT::ShadowMnTestCase<Yasp3, Dune::GDT::RealSphericalHarmonicsMomentBasis<double, double, 2, 3>,
+// false>, Dune::GDT::ShadowMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 1, 1, 3>, true>
+# endif
+// Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, false>
+// Our shifted qr eigensolver fails for this problem, needs better shifting strategy
+# if HAVE_MKL || HAVE_LAPACKE || HAVE_EIGEN
+// ,
+// Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::HatFunctionMomentBasis<double, 3, double, 0, 1, 3>, true>
+# endif
+# if HAVE_QHULL
+// ,
+// Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::PartialMomentBasis<double, 3, double, 0, 1, 3>, false>,
+// Dune::GDT::PointSourceMnTestCase<Yasp3, Dune::GDT::PartialMomentBasis<double, 3, double, 0, 1, 3>, true>
+# endif
+#endif
+ >;
+
+TYPED_TEST_CASE(HyperbolicEntropicCoordsMnTest, YaspGridTestCasesAll);
+TYPED_TEST(HyperbolicEntropicCoordsMnTest, check)
+{
+ this->run();
+}
diff --git a/dune/gdt/tools/timestepper/adaptive-rungekutta.hh b/dune/gdt/tools/timestepper/adaptive-rungekutta.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4e46bf4b1cd26eec5aee0a22ce1fb1bfec2d0628
--- /dev/null
+++ b/dune/gdt/tools/timestepper/adaptive-rungekutta.hh
@@ -0,0 +1,410 @@
+// 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_TIMESTEPPER_ADAPTIVE_RUNGEKUTTA_HH
+#define DUNE_GDT_TIMESTEPPER_ADAPTIVE_RUNGEKUTTA_HH
+
+#include <utility>
+
+#include <dune/gdt/operators/interfaces.hh>
+
+#include <dune/xt/common/memory.hh>
+#include <dune/xt/common/string.hh>
+
+#include "interface.hh"
+
+
+namespace Dune {
+namespace GDT {
+
+
+namespace internal {
+
+
+// unspecialized
+template <class RangeFieldType, TimeStepperMethods method>
+struct AdaptiveButcherArrayProvider
+{
+ static_assert(AlwaysFalse<RangeFieldType>::value,
+ "You cannot use AdaptiveRungeKuttaTimeStepper with this value of TimeStepperMethods!");
+};
+
+// user-provided Butcher array
+template <class RangeFieldType>
+struct AdaptiveButcherArrayProvider<RangeFieldType, TimeStepperMethods::adaptive_rungekutta_other>
+{
+ static Dune::DynamicMatrix<RangeFieldType> A()
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "You have to provide a Butcher array in AdaptiveRungeKuttaTimeStepper's constructor for this method!");
+ return Dune::DynamicMatrix<RangeFieldType>();
+ }
+
+ static Dune::DynamicVector<RangeFieldType> b_1()
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "You have to provide a Butcher array in AdaptiveRungeKuttaTimeStepper's constructor for this method!");
+ return Dune::DynamicVector<RangeFieldType>();
+ }
+
+ static Dune::DynamicVector<RangeFieldType> b_2()
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "You have to provide a Butcher array in AdaptiveRungeKuttaTimeStepper's constructor for this method!");
+ return Dune::DynamicVector<RangeFieldType>();
+ }
+
+ static Dune::DynamicVector<RangeFieldType> c()
+ {
+ DUNE_THROW(Dune::NotImplemented,
+ "You have to provide a Butcher array in AdaptiveRungeKuttaTimeStepper's constructor for this method!");
+ return Dune::DynamicVector<RangeFieldType>();
+ }
+};
+
+// Bogacki-Shampine (adaptive RK23)
+template <class RangeFieldType>
+class AdaptiveButcherArrayProvider<RangeFieldType, TimeStepperMethods::bogacki_shampine>
+{
+public:
+ static Dune::DynamicMatrix<RangeFieldType> A()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicMatrix<RangeFieldType>>(
+ "[0 0 0 0; 0.5 0 0 0; 0 0.75 0 0; " + Dune::XT::Common::to_string(2.0 / 9.0, 15) + " "
+ + Dune::XT::Common::to_string(1.0 / 3.0, 15) + " " + Dune::XT::Common::to_string(4.0 / 9.0, 15) + " 0]");
+ }
+
+ static Dune::DynamicVector<RangeFieldType> b_1()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicVector<RangeFieldType>>(
+ "[" + Dune::XT::Common::to_string(2.0 / 9.0, 15) + " " + Dune::XT::Common::to_string(1.0 / 3.0, 15) + " "
+ + Dune::XT::Common::to_string(4.0 / 9.0, 15) + " 0]");
+ }
+
+ static Dune::DynamicVector<RangeFieldType> b_2()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicVector<RangeFieldType>>(
+ "[" + Dune::XT::Common::to_string(7.0 / 24.0, 15) + " " + Dune::XT::Common::to_string(1.0 / 4.0, 15) + " "
+ + Dune::XT::Common::to_string(1.0 / 3.0, 15) + " " + Dune::XT::Common::to_string(1.0 / 8.0, 15) + " 0]");
+ }
+
+ static Dune::DynamicVector<RangeFieldType> c()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicVector<RangeFieldType>>("[0.5 0.75 1 0]");
+ }
+};
+
+// Dormand-Prince (adaptive RK45)
+template <class RangeFieldType>
+class AdaptiveButcherArrayProvider<RangeFieldType, TimeStepperMethods::dormand_prince>
+{
+public:
+ static Dune::DynamicMatrix<RangeFieldType> A()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicMatrix<RangeFieldType>>(
+ std::string("[0 0 0 0 0 0 0;") + " 0.2 0 0 0 0 0 0;" + " 0.075 0.225 0 0 0 0 0;" + " "
+ + Dune::XT::Common::to_string(44.0 / 45.0, 15) + " " + Dune::XT::Common::to_string(-56.0 / 15.0, 15) + " "
+ + Dune::XT::Common::to_string(32.0 / 9.0, 15) + " 0 0 0 0;" + " "
+ + Dune::XT::Common::to_string(19372.0 / 6561.0, 15) + " " + Dune::XT::Common::to_string(-25360.0 / 2187.0, 15)
+ + " " + Dune::XT::Common::to_string(64448.0 / 6561.0, 15) + " "
+ + Dune::XT::Common::to_string(-212.0 / 729.0, 15) + " 0 0 0;" + " "
+ + Dune::XT::Common::to_string(9017.0 / 3168.0, 15) + " " + Dune::XT::Common::to_string(-355.0 / 33.0, 15) + " "
+ + Dune::XT::Common::to_string(46732.0 / 5247.0, 15) + " " + Dune::XT::Common::to_string(49.0 / 176.0, 15) + " "
+ + Dune::XT::Common::to_string(-5103.0 / 18656.0, 15) + " 0 0;" + " "
+ + Dune::XT::Common::to_string(35.0 / 384.0, 15) + " 0 " + Dune::XT::Common::to_string(500.0 / 1113.0, 15) + " "
+ + Dune::XT::Common::to_string(125.0 / 192.0, 15) + " " + Dune::XT::Common::to_string(-2187.0 / 6784.0, 15) + " "
+ + Dune::XT::Common::to_string(11.0 / 84.0, 15) + " 0]");
+ }
+
+ static Dune::DynamicVector<RangeFieldType> b_1()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicVector<RangeFieldType>>(
+ "[" + Dune::XT::Common::to_string(35.0 / 384.0, 15) + " 0 " + Dune::XT::Common::to_string(500.0 / 1113.0, 15)
+ + " " + Dune::XT::Common::to_string(125.0 / 192.0, 15) + " " + Dune::XT::Common::to_string(-2187.0 / 6784.0, 15)
+ + " " + Dune::XT::Common::to_string(11.0 / 84.0, 15) + " 0]");
+ }
+
+ static Dune::DynamicVector<RangeFieldType> b_2()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicVector<RangeFieldType>>(
+ "[" + Dune::XT::Common::to_string(5179.0 / 57600.0, 15) + " 0 "
+ + Dune::XT::Common::to_string(7571.0 / 16695.0, 15) + " " + Dune::XT::Common::to_string(393.0 / 640.0, 15) + " "
+ + Dune::XT::Common::to_string(-92097.0 / 339200.0, 15) + " " + Dune::XT::Common::to_string(187.0 / 2100.0, 15)
+ + " " + Dune::XT::Common::to_string(1.0 / 40.0, 15) + "]");
+ }
+
+ static Dune::DynamicVector<RangeFieldType> c()
+ {
+ return Dune::XT::Common::from_string<Dune::DynamicVector<RangeFieldType>>(
+ "[0 0.2 0.3 0.8 " + Dune::XT::Common::to_string(8.0 / 9.0, 15) + " 1 1]");
+ }
+}; // Dormand-Prince (RK45)
+
+
+} // namespace internal
+
+
+/** \brief Time stepper using adaptive Runge Kutta methods
+ *
+ * Timestepper using adaptive Runge Kutta methods to solve equations of the form u_t = r * L(u, t) where u is a
+ * discrete function, L an operator acting on u and \alpha a scalar factor (e.g. -1).
+ * The specific Runge Kutta method can be chosen as the third template argument. If your desired Runge Kutta method is
+ * not contained in AdaptiveRungeKuttaMethods, choose AdaptiveRungeKuttaMethods::other and
+ * supply a DynamicMatrix< RangeFieldType > A and vectors b_1, b_2 (DynamicVector< RangeFieldType >) and c
+ * (DynamicVector< RangeFieldType >) in the constructor. Here, A, b_1, b_2 and c form the butcher tableau (see
+ * https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods#Embedded_methods, A is composed of the coefficients
+ * a_{ij}, b_1 of b_j, b_2 of b_j^* and c of c_j). The default is the Dormand-Prince RK45 method.
+ * In each time step, the error is estimated using the difference between the two solutions obtained using either b_1 or
+ * b_2. If the estimated error is higher than a specified tolerance tol, the calculation is repeated with a smaller
+ * time step. The tolerance tol and the error estimate are also used to estimate the optimal time step length for the
+ * next time step via dt_new = dt_old*min(max(0.9*(tol/error)^(1/5), scale_factor_min), scale_factor_max_);
+ *
+ * Notation: For an s-stage method,
+ * \mathbf{u}^{n+1} = \mathbf{u}^n + dt \sum_{i=0}^{s-1} b_i \mathbf{k}_i
+ * \mathbf{k}_i = L(\mathbf{u}_i, t^n + dt c_i)
+ * \mathbf{u}_i = \mathbf{u}^n + dt \sum_{j=0}^{i-1} a_{ij} \mathbf{k}_j,
+ *
+ * \tparam OperatorImp Type of operator L
+ * \tparam DiscreteFunctionImp Type of initial values and solution at a fixed time
+ * \tparam method Adaptive Runge-Kutta method that is used (default is AdaptiveRungeKuttaMethods::dormand_prince)
+ */
+template <class OperatorImp, class DiscreteFunctionImp, TimeStepperMethods method = TimeStepperMethods::dormand_prince>
+class AdaptiveRungeKuttaTimeStepper : public TimeStepperInterface<DiscreteFunctionImp>
+{
+ typedef TimeStepperInterface<DiscreteFunctionImp> BaseType;
+ typedef typename internal::AdaptiveButcherArrayProvider<typename BaseType::RangeFieldType, method>
+ ButcherArrayProviderType;
+
+public:
+ typedef OperatorImp OperatorType;
+ typedef DiscreteFunctionImp DiscreteFunctionType;
+
+ typedef typename DiscreteFunctionType::DomainFieldType DomainFieldType;
+ typedef typename DiscreteFunctionType::RangeFieldType RangeFieldType;
+ typedef typename Dune::DynamicMatrix<RangeFieldType> MatrixType;
+ typedef typename Dune::DynamicVector<RangeFieldType> VectorType;
+ typedef typename std::vector<std::pair<RangeFieldType, DiscreteFunctionType>> SolutionType;
+
+ /**
+ * \brief Constructor for AdaptiveRungeKuttaTimeStepper time stepper
+ * \param op Operator L
+ * \param initial_values Discrete function containing initial values for u at time t_0.
+ * \param r Scalar factor (see above, default is 1)
+ * \param t_0 Initial time (default is 0)
+ * \param tol Error tolerance for the adaptive scheme (default is 1e-4)
+ * \param scale_factor_min Minimum allowed factor for time step scaling (default is 0.2).
+ * \param scale_factor_max Maximum allowed factor for time step scaling (default is 5).
+ * \param A Coefficient matrix (only provide if you use AdaptiveRungeKuttaMethods::other)
+ * \param b_1 First set of coefficients (only provide if you use AdaptiveRungeKuttaMethods::other)
+ * \param b_2 Second set of coefficients (only provide if you use AdaptiveRungeKuttaMethods::other)
+ * \param c Coefficients for time steps (only provide if you use AdaptiveRungeKuttaMethods::other)
+ */
+ AdaptiveRungeKuttaTimeStepper(const OperatorType& op,
+ DiscreteFunctionType& initial_values,
+ const RangeFieldType r = 1.0,
+ const double t_0 = 0.0,
+ const RangeFieldType tol = 1e-4,
+ const RangeFieldType scale_factor_min = 0.2,
+ const RangeFieldType scale_factor_max = 5,
+ const MatrixType& A = ButcherArrayProviderType::A(),
+ const VectorType& b_1 = ButcherArrayProviderType::b_1(),
+ const VectorType& b_2 = ButcherArrayProviderType::b_2(),
+ const VectorType& c = ButcherArrayProviderType::c())
+ : BaseType(t_0, initial_values)
+ , op_(op)
+ , r_(r)
+ , tol_(tol)
+ , scale_factor_min_(scale_factor_min)
+ , scale_factor_max_(scale_factor_max)
+ , u_tmp_(BaseType::current_solution())
+ , A_(A)
+ , b_1_(b_1)
+ , b_2_(b_2)
+ , c_(c)
+ , b_diff_(b_2_ - b_1_)
+ , num_stages_(A_.rows())
+ {
+ assert(Dune::XT::Common::FloatCmp::gt(tol_, 0.0));
+ assert(Dune::XT::Common::FloatCmp::le(scale_factor_min_, 1.0));
+ assert(Dune::XT::Common::FloatCmp::ge(scale_factor_max_, 1.0));
+ assert(A_.rows() == A_.cols() && "A has to be a square matrix");
+ assert(b_1_.size() == A_.rows());
+ assert(b_2_.size() == A_.rows());
+ assert(c_.size() == A_.rows());
+#ifndef NDEBUG
+ for (size_t ii = 0; ii < A_.rows(); ++ii) {
+ for (size_t jj = ii; jj < A_.cols(); ++jj) {
+ assert(Dune::XT::Common::FloatCmp::eq(A_[ii][jj], 0.0)
+ && "A has to be a lower triangular matrix with 0 on the main diagonal");
+ }
+ }
+#endif // NDEBUG
+ // store as many discrete functions as needed for intermediate stages
+ for (size_t ii = 0; ii < num_stages_; ++ii) {
+ stages_k_.emplace_back(current_solution());
+ }
+ } // constructor AdaptiveRungeKuttaTimeStepper
+
+ using BaseType::current_solution;
+ using BaseType::current_time;
+ using BaseType::solve;
+
+ virtual RangeFieldType solve(const RangeFieldType t_end,
+ const RangeFieldType initial_dt,
+ const size_t num_save_steps,
+ const size_t num_output_steps,
+ const bool save_solution,
+ const bool visualize,
+ const bool with_half_steps,
+ const bool write_discrete,
+ const bool write_exact,
+ const std::string prefix,
+ typename BaseType::DiscreteSolutionType& sol,
+ const typename BaseType::VisualizerType& visualizer,
+ const typename BaseType::StringifierType& stringifier,
+ const typename BaseType::GridFunctionType& exact_solution) override final
+ {
+ const auto ret = BaseType::solve(t_end,
+ initial_dt,
+ num_save_steps,
+ num_output_steps,
+ save_solution,
+ visualize,
+ with_half_steps,
+ write_discrete,
+ write_exact,
+ prefix,
+ sol,
+ visualizer,
+ stringifier,
+ exact_solution);
+ // in a fractional step scheme, we cannot use last_stage_of_previous_step
+ last_stage_of_previous_step_ = nullptr;
+ return ret;
+ }
+
+ RangeFieldType step(const RangeFieldType dt, const RangeFieldType max_dt) override final
+ {
+ RangeFieldType actual_dt = std::min(dt, max_dt);
+ RangeFieldType mixed_error = std::numeric_limits<RangeFieldType>::max();
+ RangeFieldType time_step_scale_factor = 1.0;
+
+ auto& t = current_time();
+ auto& u_n = current_solution();
+
+ while (Dune::XT::Common::FloatCmp::gt(mixed_error, tol_)) {
+ bool skip_error_computation = false;
+ actual_dt *= time_step_scale_factor;
+ size_t first_stage_to_compute = 0;
+ if (last_stage_of_previous_step_) {
+ stages_k_[0].dofs().vector() = last_stage_of_previous_step_->dofs().vector();
+ first_stage_to_compute = 1;
+ }
+
+ for (size_t ii = first_stage_to_compute; ii < num_stages_; ++ii) {
+ std::fill(stages_k_[ii].dofs().vector().begin(), stages_k_[ii].dofs().vector().end(), RangeFieldType(0.));
+ u_tmp_.dofs().vector() = u_n.dofs().vector();
+ for (size_t jj = 0; jj < ii; ++jj)
+ u_tmp_.dofs().vector() += stages_k_[jj].dofs().vector() * (actual_dt * r_ * (A_[ii][jj]));
+ try {
+ op_.apply(u_tmp_.dofs().vector(), stages_k_[ii].dofs().vector(), t + actual_dt * c_[ii]);
+ } catch (const Dune::MathError& e) {
+ std::cout << "Caught error " << e.what() << std::endl;
+ mixed_error = 1e10;
+ skip_error_computation = true;
+ time_step_scale_factor = 0.5;
+ break;
+ }
+ }
+
+ if (!skip_error_computation) {
+ // compute error vector
+ u_tmp_.dofs().vector() = stages_k_[0].dofs().vector() * b_diff_[0];
+ for (size_t ii = 1; ii < num_stages_; ++ii)
+ u_tmp_.dofs().vector() += stages_k_[ii].dofs().vector() * b_diff_[ii];
+ u_tmp_.dofs().vector() *= actual_dt * r_;
+
+ // calculate u at timestep n+1
+ for (size_t ii = 0; ii < num_stages_; ++ii)
+ u_n.dofs().vector() += stages_k_[ii].dofs().vector() * (actual_dt * r_ * b_1_[ii]);
+
+ // scale error, use absolute error if norm is less than 0.01 and relative error else
+ auto& diff_vector = u_tmp_.dofs().vector();
+ for (size_t ii = 0; ii < diff_vector.size(); ++ii) {
+ if (std::abs(u_n.dofs().vector()[ii]) > 0.01)
+ diff_vector[ii] /= std::abs(u_n.dofs().vector()[ii]);
+ }
+ mixed_error = diff_vector.sup_norm();
+ // scale dt to get the estimated optimal time step length, TODO: adapt formula
+ time_step_scale_factor =
+ std::min(std::max(0.9 * std::pow(tol_ / mixed_error, 1.0 / 5.0), scale_factor_min_), scale_factor_max_);
+
+ if (mixed_error > tol_) { // go back from u at timestep n+1 to timestep n
+ for (size_t ii = 0; ii < num_stages_; ++ii)
+ u_n.dofs().vector() += stages_k_[ii].dofs().vector() * (-1.0 * r_ * actual_dt * b_1_[ii]);
+ }
+ }
+ } // while (mixed_error > tol_)
+ if (!last_stage_of_previous_step_)
+ last_stage_of_previous_step_ = Dune::XT::Common::make_unique<DiscreteFunctionType>(u_n);
+ last_stage_of_previous_step_->dofs().vector() = stages_k_[num_stages_ - 1].dofs().vector();
+
+#if 0
+ const auto u_local_func = u_n.local_discrete_function();
+ for (auto&& element : Dune::elements(u_n.space().grid_view())) {
+ u_local_func->bind(element);
+ for (size_t ii = 0; ii < BaseType::dimRange; ++ii) {
+// constexpr double min_val = -100.;
+ constexpr double max_val = 1000.;
+ const auto entry_ii = u_local_func->dofs().get_entry(ii);
+ if (std::abs(entry_ii) > max_val) {
+// std::cout << "limited from " << u_local_func->dofs().get_entry(ii) << " to " << min_val << " in entity " << element.geometry().center() << std::endl;
+// std::cout << "Entries are: ";
+// for (size_t jj = 0; jj < BaseType::dimRange; ++jj) {
+// std::cout << u_local_func->dofs().get_entry(jj) << " ";
+// }
+// std::cout << std::endl;
+ last_stage_of_previous_step_ = nullptr;
+// u_local_func->dofs().set_entry(ii, min_val);
+ std::cout << "Replacing " << entry_ii << " by " << std::copysign(max_val, entry_ii) << std::endl;
+ u_local_func->dofs().set_entry(ii, std::copysign(max_val, entry_ii));
+ }
+ } // ii
+ } // elements
+#endif
+
+ t += actual_dt;
+
+ return actual_dt * time_step_scale_factor;
+ } // ... step(...)
+
+private:
+ const OperatorType& op_;
+ const RangeFieldType r_;
+ const RangeFieldType tol_;
+ const RangeFieldType scale_factor_min_;
+ const RangeFieldType scale_factor_max_;
+ DiscreteFunctionType u_tmp_;
+ const MatrixType A_;
+ const VectorType b_1_;
+ const VectorType b_2_;
+ const VectorType c_;
+ const VectorType b_diff_;
+ std::vector<DiscreteFunctionType> stages_k_;
+ const size_t num_stages_;
+ std::unique_ptr<DiscreteFunctionType> last_stage_of_previous_step_;
+}; // class AdaptiveRungeKuttaTimeStepper
+
+
+} // namespace GDT
+} // namespace Dune
+
+#endif // DUNE_GDT_TIMESTEPPER_ADAPTIVE_RUNGEKUTTA_HH